Code
library(here)
library(flextable)
library(plotly)
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("02_scripts/03b_fun_calib_ntd.R"))
source(here("02_scripts/06_interventions_models.R"))Final Project | EPIB676
library(here)
library(flextable)
library(plotly)
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("02_scripts/03b_fun_calib_ntd.R"))
source(here("02_scripts/06_interventions_models.R"))—What is the decision problem your model aims to inform? What is already known, and what are the key questions you hope to address? Why is this important?—–
Opioid overdose deaths and related harms epidemic began in Canada in 2016 and is still ongoing. As a result, through the Substance Use and Addictions Program (SUAP) and other related programs, the government is funding projects that support people who use drugs and aim to mitigate the harms and deaths resulting from opioid overdose epidemic. For this project, I aim to evaluate the future economic impacts of two of these projects and their effect on reducing opioid-related deaths and harms.
This project aims to evaluate the future economic impact of these interventions by conducting a cost-benefit analysis to examine the difference in healthcare costs between the interventions and their projected effect on reducing opioid-related deaths and other important health outcomes
This project is a team project, I worked on the modeling part but many of the decisions and the work upon which the modeling was built was done by other teams members (DP, MS, IW). The model was designed and conceptualized by DP. A systematic review was conducted by DP, MS and IW to extract the effects of the interventions from the literature. Transition probabilities, and healthcare costs were extracted by DP from either the literature or were estimated based on scientific assumptions in consultation with experts working in the field and/or people with lived experiences. Many of the decisions made for the modelling part were made in consultation with DP.
An open cohort Markov model was used (fig. 1) that consists of 31 states representing the pathways of opioid use in Canada over 15 years (2015 to 2029) with monthly cycles for 15+ Canadian population (designed and conceptualized by DP). It is an open cohort model that takes into account the changes in population over time by adding people to the first state (Pain free, no use) every cycle. The yearly increase was assumed to be constant throughout the year (meaning that each increase was divided by 12, as for years in the future, projections of population size in Canada were used and proportion of those 15+ was estimated based on previous data, then a similar approach was used to estimate the monthly increase in the cohort).
Table 1: List of variables corresponding to the states in the model
flextable(ini_states_tbl_org %>%
select(`Description `, `Variable `)) |>
autofit()Description | Variable |
|---|---|
Pain free, no use | BN_PN |
Acute pain, no use | BN_ACUTE |
Chronic (non-cancer) pain, no use | BN_CHRONIC |
Cancer, no use | BN_CANCER |
Other (i.e., cough, diarrhea, sickle cell disease), no use | BN_OTHER |
Palliative, no use | BN_PALLIATIVE |
Acute pain, Rx use | BPO_ACUTE |
Chronic (non-cancer) pain, Rx use | BPO_CHRONIC |
Cancer, Rx use | BPO_CANCER |
Other (i.e., cough, diarrhea, sickle cell pain), Rx use | BPO_OTHER |
Palliative, Rx use | BPO_PALLIATIVE |
Rx opioid misuse | BPO_MISUSE |
Illicit opioid use | BI_ILLICIT |
Detox / Withdrawal Management | BS_DETOX |
OAT initiation | BS_OAT_INI |
OAT maintenance | BS_OAT_MAINT |
OAT / Safe Supply | BS_OAT_SS |
Safe Supply | BS_SS |
Overdoes - illicit | BO_OD_ILLICIT |
Overdose - misuse | BO_OD_RX |
Moderate brain injury | BO_MOD_BI |
Severe brain injury | BO_SEVERE_BI |
Severe brain injury Out | BO_SEVERE_BI_OUT |
Opioid-related death | BO_OD_DEATH |
Death | BO_DEATH |
R: Illicit opioid use | BR_ILLICIT |
R: OAT initiation | BR_OAT_INI |
R: OAT maintenance | BR_OAT_MAINT |
R: OAT / Safe Supply | BR_OAT_SS |
R: Safe Supply | BR_SS |
R: Illicit opioid overdose | BR_OD_ILLICIT |
Characteristics of this model:
This model was run for 5 scenarios in regards to the previously mentioned interventions: 1) No interventions, 2) Naloxone only, 3) Prescription guidelines only, 4) Safer Supply only, and 5) all 3 interventions. The effects of each of the interventions were estimated by other team members who conducted a systematic review of the literature (DP, MS, IW). Transition probabilities were changed for the naloxone and prescription guidelines interventions (as a result the 1-rest transition probability will be changed).
For the naloxone intervention the following transition probabilities were changed starting from January 2017:
For the prescription guidelines the following transition probabilities were changed starting from May 2017:
For safer supply to account for the effect of the pilot programs introduced in 2016, 1000 people were moved from illicit drug use state to OAT/Safer supply state in January 2016. Then to account for compassionate prescribing due to the COVID-19 panademic 10,000 people were moved from illicit drug use state to OAT/Safer supply state in April 2020. Before January 2016 all transitioning into and out of safer supply states (Safer Supply, OAT/Safer Supply, Repeat OAT/Safer Supply and Repeat Safer Supply). Starting from January 2016, all transitioning out of OAT/Safer Supply (and staying in) were changed to their corresponding non-zero values. Until in 2023, when the intervention is supposed to be implemented, all transitioning probabilities into and out of safer supply states were changed to their non-zero values.
For the cost-benefit analysis the effects of these interventions individually and all together was on interest, therefore 4 different comparisons were made: 1) Naloxone vs No Interventions, 2) Safer Supply vs No Interventions, 3) Prescription Guidelines vs No Interventions and 4) All Interventions vs No Interventions.
Primary Outcomes:
Secondary Outcomes - more epidemiological measures:
For the primary outcomes changes and percent changes were calculated for the 4 different comparison groups while the secondary outcomes were calculated for each of th 5 scenarios.
After designing the model structure, initial populations, and transitions, we calibrated the model using data from 2015 to 2021. Our calibration targets were total deaths (excluding opioid-related overdose deaths), the total number of people in Opioid Agonist Treatment (OAT), and proportions of prescription opioid use. To obtain this data, we used various sources, including Statistics Canada, the Canadian Institute for Health Information’s report on Opioid prescribing in Canada, and scientific literature (DP).
Even though we had opioid-related OD deaths numbers, we didn’t use them as targets because they are directly affected by different interventions that has been implemented so we couldn’t calibrate the no intervention parameters to them. However, we used them to assess if the estimates of opioid-related OD deaths from the model of the no intervention scenario are sound (they should be more than the targets).
We selected parameters to calibrate based on a one-way sensitivity analysis, where we changed one transition probability at a time, using lower and upper ranges, and observed how this impacted predicted outcomes (overall deaths and OD deaths). The one-way sensitivity analysis was conducted twice, 1) accounting for the fentanyl contamination (Scenario 1) and 2) without accounting for fentanyl contamination (Scenario 2). We included transition probabilities that resulted in a notable change in predicted outcomes (defined as 2.5%) in the set of parameters to be calibrated, in addition to all transition probabilities to death or opioid-related overdose deaths (a total of 38 parameters for scenario 1 and 39 paramters for scenario 2). We then conducted a Maximum A Posteriori (MAP) estimation using Nelder Mead, a sampling method that aims to find the global optimum twice once with the model for both scenarios. This allowed us to select the parameter set that best approximated the targets and provided us with information about changes in the prior information for some of the parameters. Overall, we updated the following baseline transition probabilities:
Then we used the Sample Importance Resampling (SIR) Bayesian calibration method after updating the prior information based on the MAP results, and we ran the model using the sample sets (that were randomly sampled from the priors defined for the parameters), estimated likelihood for each set, and weighted the sets by likelihood to approximate posteriors.
Probabilistic sensitivity analysis (PSA) was conducted on primary outcomes to account for the uncertainty in our result. The SIR generated posteriors for the parameters that were calibrated were used as priors and for other transition probabilities, a uniform distribution was used with ranges either from the literature and if not available, +/-5% were used as ranges. From these priors probabilistic samples were randomly samples and the analysis was simulated 10 000 times.
What are your main findings? Should estimate some outcome(s) for >1 alternative, and potentially conduct incremental analysis (compute ICERS) depending on your framework. Report uncertainty analysis findings to give some idea of the robustness of your findings
The model approximated the calibration targets well: even though the distribution of prevalence of opioid prescription did not include the target values, the difference between the values was negligible. The distribution for overall deaths included the target values and the distribution for opioid-related deaths was larger than the target values (what we wanted the model to result in). However the distribution of OAT counts doesn’t include the target values and that might be either because these values aren’t accurate enough to represent the values estimated by the model or that the calibration process didn’t work as well. Overall, the estimates from the model seem to be approximating values found in the literature.
m_outcomes_psa <- readRDS(here("01_data/m_outcomes_psa.RDS"))
point_est <- readRDS(here("01_data/point_estiamte.RDS"))
total_death_oddeaths <- readRDS(here("01_data/total_deaths_oddeaths.RDS"))
mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)
mean_cost_diff_nalox <- mean(m_outcomes_psa[, "cost_diff_nalox"])
ci_cost_diff_nalox <- quantile(m_outcomes_psa[, "cost_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_ss <- mean(m_outcomes_psa[, "cost_diff_ss"])
ci_cost_diff_ss <- quantile(m_outcomes_psa[, "cost_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_pg <- mean(m_outcomes_psa[, "cost_diff_pg"])
ci_cost_diff_pg <- quantile(m_outcomes_psa[, "cost_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_all <- mean(m_outcomes_psa[, "cost_diff_all"])
ci_cost_diff_all <- quantile(m_outcomes_psa[, "cost_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff <- c(paste0(round(ci_cost_diff_nalox[3]/1000000, 0), " [",
round(ci_cost_diff_nalox[1]/1000000, 0),", ",
round(ci_cost_diff_nalox[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_ss[3]/1000000, 0), " [",
round(ci_cost_diff_ss[1]/1000000, 0),", ",
round(ci_cost_diff_ss[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_pg[3]/1000000, 0), " [",
round(ci_cost_diff_pg[1]/1000000, 0),", ",
round(ci_cost_diff_pg[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_all[3]/1000000, 0), " [",
round(ci_cost_diff_all[1]/1000000, 0),", ",
round(ci_cost_diff_all[2]/1000000, 0), "]"))
mean_death_diff_nalox <- mean(m_outcomes_psa[, "death_diff_nalox"])
ci_death_diff_nalox <- quantile(m_outcomes_psa[, "death_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_ss <- mean(m_outcomes_psa[, "death_diff_ss"])
ci_death_diff_ss <- quantile(m_outcomes_psa[, "death_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_pg <- mean(m_outcomes_psa[, "death_diff_pg"])
ci_death_diff_pg <- quantile(m_outcomes_psa[, "death_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_all <- mean(m_outcomes_psa[, "death_diff_all"])
ci_death_diff_all <- quantile(m_outcomes_psa[, "death_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff <- c(paste0(round(ci_death_diff_nalox[3], 0), " [",
round(ci_death_diff_nalox[1], 0),", ",
round(ci_death_diff_nalox[2], 0), "]"),
paste0(round(ci_death_diff_ss[3], 0), " [",
round(ci_death_diff_ss[1], 0),", ",
round(ci_death_diff_ss[2], 0), "]"),
paste0(round(ci_death_diff_pg[3], 0), " [",
round(ci_death_diff_pg[1], 0),", ",
round(ci_death_diff_pg[2], 0), "]"),
paste0(round(ci_death_diff_all[3], 0), " [",
round(ci_death_diff_all[1], 0),", ",
round(ci_death_diff_all[2], 0), "]"))
mean_oddeath_diff_nalox <- mean(m_outcomes_psa[, "oddeath_diff_nalox"])
ci_oddeath_diff_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_ss <- mean(m_outcomes_psa[, "oddeath_diff_ss"])
ci_oddeath_diff_ss <- quantile(m_outcomes_psa[, "oddeath_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_pg <- mean(m_outcomes_psa[, "oddeath_diff_pg"])
ci_oddeath_diff_pg <- quantile(m_outcomes_psa[, "oddeath_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_all <- mean(m_outcomes_psa[, "oddeath_diff_all"])
ci_oddeath_diff_all <- quantile(m_outcomes_psa[, "oddeath_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff <- c(paste0(round(ci_oddeath_diff_nalox[3], 0), " [",
round(ci_oddeath_diff_nalox[1], 0),", ",
round(ci_oddeath_diff_nalox[2], 0), "]"),
paste0(round(ci_oddeath_diff_ss[3], 0), " [",
round(ci_oddeath_diff_ss[1], 0),", ",
round(ci_oddeath_diff_ss[2], 0), "]"),
paste0(round(ci_oddeath_diff_pg[3], 0), " [",
round(ci_oddeath_diff_pg[1], 0),", ",
round(ci_oddeath_diff_pg[2], 0), "]"),
paste0(round(ci_oddeath_diff_all[3], 0), " [",
round(ci_oddeath_diff_all[1], 0),", ",
round(ci_oddeath_diff_all[2], 0), "]"))
effects_tbl <- cbind.data.frame(mod_names[-1], mean_cost_diff,
mean_death_diff, mean_oddeath_diff)
effects_tbl |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
mean_cost_diff = "Discounted Net Present \n Costs in Millions ",
mean_death_diff = "Deaths ",
mean_oddeath_diff = "Opioid-related Overdose Deaths")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present | Deaths | Opioid-related Overdose Deaths |
Naloxone | 209 [107, 362] | 625 [414, 926] | -10089 [-15009, -6697] |
Safer Supply | 4216 [3651, 4821] | -14180 [-17532, -11074] | -3221 [-4105, -2529] |
Prescription Guidelines | -25455 [-31480, -19866] | 14417 [8610, 20565] | -5458 [-11045, -1624] |
All Interventions | -21058 [-27059, -15403] | 960 [-5687, 7753] | -18544 [-28496, -11309] |
mean_cost_diff_per_nalox <- mean(m_outcomes_psa[, "cost_diff_per_nalox"])
ci_cost_diff_per_nalox <- quantile(m_outcomes_psa[, "cost_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_ss <- mean(m_outcomes_psa[, "cost_diff_per_ss"])
ci_cost_diff_per_ss <- quantile(m_outcomes_psa[, "cost_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_pg <- mean(m_outcomes_psa[, "cost_diff_per_pg"])
ci_cost_diff_per_pg <- quantile(m_outcomes_psa[, "cost_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_all <- mean(m_outcomes_psa[, "cost_diff_per_all"])
ci_cost_diff_per_all <- quantile(m_outcomes_psa[, "cost_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per <- c(paste0(round(ci_cost_diff_per_nalox[3], 2), " [",
round(ci_cost_diff_per_nalox[1], 2),", ",
round(ci_cost_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_cost_diff_per_ss[3], 2), " [",
round(ci_cost_diff_per_ss[1], 2),", ",
round(ci_cost_diff_per_ss[2], 2), "]%"),
paste0(round(ci_cost_diff_per_pg[3], 2), " [",
round(ci_cost_diff_per_pg[1], 2),", ",
round(ci_cost_diff_per_pg[2], 2), "]%"),
paste0(round(ci_cost_diff_per_all[3], 2), " [",
round(ci_cost_diff_per_all[1], 2),", ",
round(ci_cost_diff_per_all[2], 2), "]%"))
mean_death_diff_per_nalox <- mean(m_outcomes_psa[, "death_diff_per_nalox"])
ci_death_diff_per_nalox <- quantile(m_outcomes_psa[, "death_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_ss <- mean(m_outcomes_psa[, "death_diff_per_ss"])
ci_death_diff_per_ss <- quantile(m_outcomes_psa[, "death_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_pg <- mean(m_outcomes_psa[, "death_diff_per_pg"])
ci_death_diff_per_pg <- quantile(m_outcomes_psa[, "death_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_all <- mean(m_outcomes_psa[, "death_diff_per_all"])
ci_death_diff_per_all <- quantile(m_outcomes_psa[, "death_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per <- c(paste0(round(ci_death_diff_per_nalox[3], 2), " [",
round(ci_death_diff_per_nalox[1], 2),", ",
round(ci_death_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_death_diff_per_ss[3], 2), " [",
round(ci_death_diff_per_ss[1], 2),", ",
round(ci_death_diff_per_ss[2], 2), "]%"),
paste0(round(ci_death_diff_per_pg[3], 2), " [",
round(ci_death_diff_per_pg[1], 2),", ",
round(ci_death_diff_per_pg[2], 2), "]%"),
paste0(round(ci_death_diff_per_all[3], 2), " [",
round(ci_death_diff_per_all[1], 2),", ",
round(ci_death_diff_per_all[2], 2), "]%"))
mean_oddeath_diff_per_nalox <- mean(m_outcomes_psa[, "oddeath_diff_per_nalox"])
ci_oddeath_diff_per_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_ss <- mean(m_outcomes_psa[, "oddeath_diff_per_ss"])
ci_oddeath_diff_per_ss <- quantile(m_outcomes_psa[, "oddeath_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_pg <- mean(m_outcomes_psa[, "oddeath_diff_per_pg"])
ci_oddeath_diff_per_pg <- quantile(m_outcomes_psa[, "oddeath_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_all <- mean(m_outcomes_psa[, "oddeath_diff_per_all"])
ci_oddeath_diff_per_all <- quantile(m_outcomes_psa[, "oddeath_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per <- c(paste0(round(ci_oddeath_diff_per_nalox[3], 2), " [",
round(ci_oddeath_diff_per_nalox[1], 2),", ",
round(ci_oddeath_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_ss[3], 2), " [",
round(ci_oddeath_diff_per_ss[1], 2),", ",
round(ci_oddeath_diff_per_ss[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_pg[3], 2), " [",
round(ci_oddeath_diff_per_pg[1], 2),", ",
round(ci_oddeath_diff_per_pg[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_all[3], 2), " [",
round(ci_oddeath_diff_per_all[1], 2),", ",
round(ci_oddeath_diff_per_all[2], 2), "]%"))
effects_tbl_per <- cbind.data.frame(mod_names[-1], mean_cost_diff_per,
mean_death_diff_per, mean_oddeath_diff_per)
effects_tbl_per |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Percent Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
mean_cost_diff_per = "Discounted Net Present Costs",
mean_death_diff_per = "Deaths ",
mean_oddeath_diff_per = "Opioid-related Overdose Deaths")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Percent Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths | Opioid-related Overdose Deaths |
Naloxone | 0.03 [0.01, 0.05]% | 0.01 [0.01, 0.02]% | -3.46 [-4.13, -3.13]% |
Safer Supply | 0.54 [0.46, 0.61]% | -0.31 [-0.38, -0.24]% | -1.11 [-1.97, -0.66]% |
Prescription Guidelines | -3.23 [-3.94, -2.55]% | 0.32 [0.19, 0.45]% | -1.81 [-2.76, -0.93]% |
All Interventions | -2.67 [-3.39, -1.98]% | 0.02 [-0.13, 0.17]% | -6.43 [-7.16, -5.58]% |
tbl_effect_plot1 <- cbind.data.frame(value = c("lb", "ub", "median"),
ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
ci_cost_diff_per_pg, ci_cost_diff_per_all,
ci_death_diff_per_nalox, ci_death_diff_per_ss,
ci_death_diff_per_pg, ci_death_diff_per_all,
ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>%
pivot_longer(2:13, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3)) %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions", Intervention)))))
tbl_effect_plot1$Intervention <- factor(tbl_effect_plot1$Intervention, levels = mod_names, labels = mod_names)
ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = median))+
geom_jitter(aes(color = Intervention, shape = Intervention), size = 2, position = position_dodge(0.25)) +
# geom_errorbar(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), width = 0.3, position = position_dodge(0.1)) +
geom_pointrange(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), position = position_dodge(0.25)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16, 17, 18)) +
# scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","Opioid-related Overdose Deaths"))+
theme_minimal())tbl_effect_plot2 <- cbind.data.frame(value = c("cost_lb", "cost_ub", "cost_median"),
ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
ci_cost_diff_per_pg, ci_cost_diff_per_all) %>%
pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3, group)) %>%
full_join(., cbind.data.frame(value = c("oddeath_lb", "oddeath_ub", "oddeath_median"),
ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>%
pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3, group)), by = "Intervention") %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions", Intervention)))))
tbl_effect_plot2$Intervention <- factor(tbl_effect_plot2$Intervention, levels = mod_names, labels = mod_names)
ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_median, y = oddeath_median)) +
geom_point(aes(color = Intervention), size = 1) +
geom_pointrange(aes(ymin = oddeath_lb, ymax = oddeath_ub, color = Intervention)) +
geom_errorbarh(aes(xmax = cost_lb, xmin = cost_ub,color = Intervention), height = 0) +
scale_color_brewer(palette = "Dark2") +
# scale_shape_manual(values = c(15, 16, 17, 18)) +
geom_hline(yintercept = 0, linewidth = 0.25) +
geom_vline(xintercept = 0, linewidth = 0.25) +
xlab("Change in Costs Compared to No Intervention (%)") +
ylab("Change in Opioid-related Overdose Deaths \n Compared to No Intervention (%)") +
theme_minimal())The results show the 3 interventions individually and all together reduced the total number of opioid-related deaths compared to no interventions, with all interventions combined resulting in -6.413539 [-5.58, -7.16] % difference compared to no intervention scenario. All interventions and prescription guidelines resulted in reduction costs as well, however Naloxone and Safer Supply increased costs slightly. Implementing all interventions reduces both costs and opioid-related OD deaths.
Regarding overall deaths though, except for safer supply, all other scenarios either increase overall deaths or result in inconclusive results.
v_yrs_prev <- c(2015:2029)
yrs_prev <- length(v_yrs_prev)
prev_fun <- function(mod, i, v_states){
prop_uncalib <- data.frame(matrix(NA,byrow = T, nrow = 12,
ncol = yrs_prev,
list(NULL,
paste0("prop_",
str_sub(v_yrs_prev, 3, 4)))))
tot_pop_uncalib_tbl <- data.frame(matrix(NA,byrow = T, nrow = 12,
ncol = yrs_prev,
list(NULL,
paste0("tot_pop_",
str_sub(v_yrs_prev, 3, 4)))))
for (j in 1:yrs_prev) {
yr <- (v_yrs_prev[1] - 1) + j
if (length(v_states) > 1)
{
prop_uncalib[, j] <- assign(
paste0("prop_", str_sub(yr, 3, 4)),
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
v_states]) /
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
} else {
prop_uncalib[, j] <- assign(
paste0("prop_", str_sub(yr, 3, 4)),
mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
v_states]) /
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
}
tot_pop_uncalib_tbl[, j] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
# final table for monthly prevalence over 15 years
prop_uncalib_tbl <- prop_uncalib %>%
pivot_longer(1:15, names_to = "grp", values_to = "value") %>%
arrange(grp) %>%
mutate(year = rep(v_yrs_prev,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:15, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop)) %>%
mutate(intervention = mod_names[i])
# table with weighted yearly prevalence
prop_uncalib_wei_mean <- prop_uncalib_tbl %>%
group_by(year) %>%
summarise(value = weighted.mean(value, tot_pop_val)) %>%
ungroup() %>%
mutate(intervention = mod_names[i])
return (list(prop_uncalib_tbl = prop_uncalib_tbl,
prop_uncalib_wei_mean = prop_uncalib_wei_mean))
}
# Prevalence of severe brain injury
noint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
nalx_prev_sevbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
ss_prev_sevbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
pg_prev_sevbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
allint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
prev_sevbi_yr_tbl <- noint_prev_sevbi_yr_tbl %>%
bind_rows(., nalx_prev_sevbi_yr_tbl, ss_prev_sevbi_yr_tbl,
pg_prev_sevbi_yr_tbl, allint_prev_sevbi_yr_tbl)Weighted Yearly Prevalence of Severe Brain Injury
ggplotly(ggplot() +
geom_point(data = prev_sevbi_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))# Prevalence of moderate Brain injury
noint_prev_modbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
nalx_prev_modbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
ss_prev_modbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
pg_prev_modbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
allint_prev_modbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
prev_modbi_yr_tbl <- noint_prev_modbi_yr_tbl %>%
bind_rows(., nalx_prev_modbi_yr_tbl, ss_prev_modbi_yr_tbl,
pg_prev_modbi_yr_tbl, allint_prev_modbi_yr_tbl)Weighted Yearly Prevalence of Moderate Brain Injury
ggplotly(ggplot() +
geom_point(data = prev_modbi_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))# Prevalence of substance use treatment
v_states_prev_tx <- c("BS_DETOX", "BS_OAT_INI", "BS_OAT_MAINT",
"BS_OAT_SS", "BS_SS", "BR_OAT_INI", "BR_OAT_MAINT",
"BR_OAT_SS", "BR_SS")
noint_prev_tx_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
nalx_prev_tx_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
ss_prev_tx_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
pg_prev_tx_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
allint_prev_tx_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
prev_tx_yr_tbl <- noint_prev_tx_yr_tbl %>%
bind_rows(., nalx_prev_tx_yr_tbl, ss_prev_tx_yr_tbl,
pg_prev_tx_yr_tbl, allint_prev_tx_yr_tbl)Weighted Yearly Prevalence of substance use treatment
ggplotly(ggplot() +
geom_point(data = prev_tx_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))# Prevalence of prescription opioid use
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target) %>%
rename(value = target) %>%
mutate(year_mon = paste(year, month.abb[7], sep = "_"))
v_states_prev_rx_opd <- c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")
noint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
nalx_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
ss_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
pg_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
allint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
prev_rx_opd_yr_tbl <- noint_prev_rx_opd_yr_tbl %>%
bind_rows(., nalx_prev_rx_opd_yr_tbl, ss_prev_rx_opd_yr_tbl,
pg_prev_rx_opd_yr_tbl, allint_prev_rx_opd_yr_tbl) %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, value) %>%
mutate(intervention = "Target"))Weighted Yearly Prevalence of prescription opioid use
ggplotly(ggplot() +
geom_point(data = prev_rx_opd_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))# Prevalence of illicit use
v_states_illicit <- c("BI_ILLICIT", "BR_ILLICIT")
noint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicit)$prop_uncalib_wei_mean
nalx_prev_illicituse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_illicit)$prop_uncalib_wei_mean
ss_prev_illicituse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_illicit)$prop_uncalib_wei_mean
pg_prev_illicituse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_illicit)$prop_uncalib_wei_mean
allint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_illicit)$prop_uncalib_wei_mean
prev_illicituse_yr_tbl <- noint_prev_illicituse_yr_tbl %>%
bind_rows(., nalx_prev_illicituse_yr_tbl, ss_prev_illicituse_yr_tbl,
pg_prev_illicituse_yr_tbl, allint_prev_illicituse_yr_tbl)Weighted Yearly Prevalence of illicit opioid use
ggplotly(ggplot() +
geom_point(data = prev_illicituse_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))# Prevalence of misuse
noint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
nalx_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
ss_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
pg_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
allint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
prev_rxmisuse_yr_tbl <- noint_prev_rxmisuse_yr_tbl %>%
bind_rows(., nalx_prev_rxmisuse_yr_tbl, ss_prev_rxmisuse_yr_tbl,
pg_prev_rxmisuse_yr_tbl, allint_prev_rxmisuse_yr_tbl)Weighted Yearly Prevalence of rx opioid misuse
ggplotly(ggplot() +
geom_point(data = prev_rxmisuse_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))# Prevalence of overdose
v_states_od <- c("BO_OD_RX", "BO_OD_ILLICIT", "BR_OD_ILLICIT")
noint_prev_od_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_od)$prop_uncalib_wei_mean
nalx_prev_od_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_od)$prop_uncalib_wei_mean
ss_prev_od_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_od)$prop_uncalib_wei_mean
pg_prev_od_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_od)$prop_uncalib_wei_mean
allint_prev_od_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_od)$prop_uncalib_wei_mean
prev_od_yr_tbl <- noint_prev_od_yr_tbl %>%
bind_rows(., nalx_prev_od_yr_tbl, ss_prev_od_yr_tbl,
pg_prev_od_yr_tbl, allint_prev_od_yr_tbl)Weighted Yearly Prevalence of opioid-related overdose
ggplotly(ggplot() +
geom_point(data = prev_od_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))The results show that none of the interventions affected prevalence of brain injury except for Safer Supply. Safer supply also resulted in an increase of prevalence in substance use treatment (which is expected since it is part of substance use treatment). Prescription Guidelines and All interventions resulted in a reduction in prevalence of prescription opioid use and prescription opioid misuse. No notable differences were observed for prevalence of opioid illicit drug use and opioid-related overdose.
If you were to develop this initial analysis into one fit for submission for peer reviewed publication, what steps would you take? Which aspects of the analysis would need to be improved? What are the key takeaways from your initial analysis (and/or, what might the takeaways be if a more complete analysis was conducted)? What is the relationship between your findings and the wider literature or policy landscape? What are the strengths and limitations of your analysis?
This analysis shows the promising effects of these interventions in reducing the direct harms of opioid crisis by reducing opioid-related OD deaths. From a healthcare prespective, when all interventions are implemented simultaneously they would reduce costs while saving lives. However, the results regarding overall deaths are not as promising, we are interested in examining the pathways that led to that increase in overall deaths.
One of the strengths of this analysis is that the model distinguished between chronic noncancer and cancer pain and included palliative care. It also included brain injuries as one of the often neglected sequalea of overdose. The uncalibrated model estimated values that were similar to target values (a strength of the model), however, it seemed that the calibration did not make a significant difference (a limitation). The calibration might not have worked due to non-identifiability problem since we have many parameters that we were trying to calibrate and few calibration targets. However that is due to the lack of data availability. This has been the main limitation of this project, there is a lack of available for the different parameters, many of which were based on scientific assumptions and knowledge of experts in the field. Therefore, as a next step, we are interested in 1) repeating the analysis for only Ontario since they have more data available so we can better calibrate the model parameters, 2) conduct a value of information analysis to identify the parameters that we need to collect data on and 3) use British Colombia data to try to further understand the pathways that led to increase in overall deaths and repeat the analysis for this province similar to Onatrio.
flextable(ini_states_tbl_org %>%
select(`Variable `, basevalue)) |>
autofit()Variable | basevalue |
|---|---|
BN_PN | 12,329,264 |
BN_ACUTE | 77,121 |
BN_CHRONIC | 6,000,000 |
BN_CANCER | 163,129 |
BN_OTHER | 6,003,590 |
BN_PALLIATIVE | 4,003 |
BPO_ACUTE | 216,300 |
BPO_CHRONIC | 3,530,940 |
BPO_CANCER | 123,666 |
BPO_OTHER | 594,870 |
BPO_PALLIATIVE | 14,158 |
BPO_MISUSE | 274,426 |
BI_ILLICIT | 290,063 |
BS_DETOX | 7,500 |
BS_OAT_INI | 2,461 |
BS_OAT_MAINT | 72,573 |
BS_OAT_SS | 0 |
BS_SS | 0 |
BO_OD_ILLICIT | 0 |
BO_OD_RX | 0 |
BO_MOD_BI | 0 |
BO_SEVERE_BI | 0 |
BO_SEVERE_BI_OUT | 0 |
BO_OD_DEATH | 0 |
BO_DEATH | 0 |
BR_ILLICIT | 0 |
BR_OAT_INI | 0 |
BR_OAT_MAINT | 0 |
BR_OAT_SS | 0 |
BR_SS | 0 |
BR_OD_ILLICIT | 0 |
flextable(pop_increase_tbl %>%
mutate(Year = as.character(Year)) %>%
rename("Increase per month" = increase_per_month)) |>
autofit()Year | Increase per month |
|---|---|
2015 | 18,637 |
2016 | 27,882 |
2017 | 32,458 |
2018 | 38,812 |
2019 | 41,247 |
2020 | 31,231 |
2021 | 19,700 |
2022 | 54,525 |
2023 | 48,524 |
2024 | 51,161 |
2025 | 50,779 |
2026 | 50,018 |
2027 | 49,005 |
2028 | 48,006 |
2029 | 47,094 |
flextable(trans_prob_tbl_new %>%
filter(basevalue != 0) %>%
rename("From - state" = var_from,
"To - state" = var_to,
Basevalue = basevalue)) |>
autofit()From - state | To - state | Basevalue |
|---|---|---|
BN_PN | BN_PN | 999.00000000 |
BN_PN | BN_ACUTE | 0.00165257 |
BN_PN | BN_CANCER | 0.00050000 |
BN_PN | BN_OTHER | 0.00069143 |
BN_PN | BN_PALLIATIVE | 0.00030770 |
BN_PN | BPO_ACUTE | 0.00110106 |
BN_PN | BPO_MISUSE | 0.00220011 |
BN_PN | BI_ILLICIT | 0.00002940 |
BN_PN | BO_DEATH | 0.00032167 |
BN_ACUTE | BN_PN | 999.00000000 |
BN_ACUTE | BN_CHRONIC | 0.60000000 |
BN_ACUTE | BPO_MISUSE | 0.00004680 |
BN_ACUTE | BI_ILLICIT | 0.00001040 |
BN_ACUTE | BO_DEATH | 0.00053201 |
BN_CHRONIC | BN_PN | 0.00010000 |
BN_CHRONIC | BN_CHRONIC | 999.00000000 |
BN_CHRONIC | BN_CANCER | 0.00042924 |
BN_CHRONIC | BN_PALLIATIVE | 0.00030770 |
BN_CHRONIC | BPO_CHRONIC | 0.01000000 |
BN_CHRONIC | BPO_MISUSE | 0.00728000 |
BN_CHRONIC | BI_ILLICIT | 0.00010000 |
BN_CHRONIC | BO_DEATH | 0.00040000 |
BN_CANCER | BN_PN | 0.00010000 |
BN_CANCER | BN_CHRONIC | 0.00003130 |
BN_CANCER | BN_CANCER | 999.00000000 |
BN_CANCER | BN_PALLIATIVE | 0.00066950 |
BN_CANCER | BPO_CANCER | 0.55000000 |
BN_CANCER | BO_DEATH | 0.01390000 |
BN_OTHER | BN_PN | 0.00053731 |
BN_OTHER | BN_OTHER | 999.00000000 |
BN_OTHER | BPO_OTHER | 0.00349866 |
BN_OTHER | BO_DEATH | 0.00060864 |
BN_PALLIATIVE | BN_PALLIATIVE | 0.18500000 |
BN_PALLIATIVE | BPO_PALLIATIVE | 0.44000000 |
BN_PALLIATIVE | BO_DEATH | 999.00000000 |
BPO_ACUTE | BN_PN | 999.00000000 |
BPO_ACUTE | BN_CHRONIC | 0.60000000 |
BPO_ACUTE | BPO_CHRONIC | 0.06200000 |
BPO_ACUTE | BPO_MISUSE | 0.00600000 |
BPO_ACUTE | BO_OD_RX | 0.00743025 |
BPO_ACUTE | BO_DEATH | 0.00053201 |
BPO_CHRONIC | BN_CHRONIC | 0.01750000 |
BPO_CHRONIC | BPO_CHRONIC | 999.00000000 |
BPO_CHRONIC | BPO_MISUSE | 0.00470310 |
BPO_CHRONIC | BO_OD_RX | 0.00160398 |
BPO_CHRONIC | BO_DEATH | 0.00040000 |
BPO_CANCER | BN_CANCER | 0.00874161 |
BPO_CANCER | BPO_CHRONIC | 0.00177217 |
BPO_CANCER | BPO_CANCER | 999.00000000 |
BPO_CANCER | BPO_PALLIATIVE | 0.00100425 |
BPO_CANCER | BPO_MISUSE | 0.00470310 |
BPO_CANCER | BO_OD_RX | 0.00097898 |
BPO_CANCER | BO_DEATH | 0.01390000 |
BPO_OTHER | BN_OTHER | 0.02825010 |
BPO_OTHER | BPO_OTHER | 999.00000000 |
BPO_OTHER | BPO_MISUSE | 0.00470310 |
BPO_OTHER | BO_OD_RX | 0.00097898 |
BPO_OTHER | BO_DEATH | 0.00040000 |
BPO_PALLIATIVE | BPO_PALLIATIVE | 999.00000000 |
BPO_PALLIATIVE | BO_DEATH | 0.38250000 |
BPO_MISUSE | BN_PN | 0.01000000 |
BPO_MISUSE | BPO_MISUSE | 999.00000000 |
BPO_MISUSE | BI_ILLICIT | 0.00426953 |
BPO_MISUSE | BS_DETOX | 0.00100000 |
BPO_MISUSE | BS_OAT_INI | 0.00100000 |
BPO_MISUSE | BO_OD_RX | 0.00097898 |
BPO_MISUSE | BO_DEATH | 0.00064333 |
BI_ILLICIT | BN_PN | 0.00001000 |
BI_ILLICIT | BN_CANCER | 0.00064386 |
BI_ILLICIT | BI_ILLICIT | 999.00000000 |
BI_ILLICIT | BS_DETOX | 0.01000000 |
BI_ILLICIT | BS_OAT_INI | 0.00700000 |
BI_ILLICIT | BO_OD_ILLICIT | 0.01320000 |
BI_ILLICIT | BO_DEATH | 0.00080801 |
BS_DETOX | BN_PN | 999.00000000 |
BS_DETOX | BI_ILLICIT | 0.22181118 |
BS_DETOX | BS_OAT_INI | 0.82000000 |
BS_DETOX | BO_OD_ILLICIT | 0.01431324 |
BS_DETOX | BO_DEATH | 0.00309346 |
BS_OAT_INI | BI_ILLICIT | 999.00000000 |
BS_OAT_INI | BS_OAT_MAINT | 0.91700000 |
BS_OAT_INI | BO_OD_ILLICIT | 0.04050000 |
BS_OAT_INI | BO_DEATH | 0.00150000 |
BS_OAT_MAINT | BN_PN | 0.00460000 |
BS_OAT_MAINT | BI_ILLICIT | 999.00000000 |
BS_OAT_MAINT | BS_OAT_MAINT | 0.71000000 |
BS_OAT_MAINT | BO_OD_ILLICIT | 0.00420000 |
BS_OAT_MAINT | BO_DEATH | 0.00038326 |
BS_OAT_SS | BS_OAT_SS | 999.00000000 |
BS_SS | BS_SS | 999.00000000 |
BO_OD_ILLICIT | BI_ILLICIT | 999.00000000 |
BO_OD_ILLICIT | BS_DETOX | 0.26800000 |
BO_OD_ILLICIT | BS_OAT_INI | 0.05871300 |
BO_OD_ILLICIT | BO_MOD_BI | 0.03000000 |
BO_OD_ILLICIT | BO_SEVERE_BI | 0.00086172 |
BO_OD_ILLICIT | BO_OD_DEATH | 0.01900000 |
BO_OD_RX | BPO_MISUSE | 999.00000000 |
BO_OD_RX | BS_DETOX | 0.26800000 |
BO_OD_RX | BS_OAT_INI | 0.05871300 |
BO_OD_RX | BO_MOD_BI | 0.03000000 |
BO_OD_RX | BO_SEVERE_BI | 0.00086172 |
BO_OD_RX | BO_OD_DEATH | 0.01380000 |
BO_MOD_BI | BO_DEATH | 0.00874161 |
BO_MOD_BI | BR_ILLICIT | 999.00000000 |
BO_MOD_BI | BR_OAT_INI | 0.15100000 |
BO_SEVERE_BI | BO_SEVERE_BI_OUT | 999.00000000 |
BO_SEVERE_BI | BO_DEATH | 0.67000000 |
BO_SEVERE_BI_OUT | BO_SEVERE_BI_OUT | 999.00000000 |
BO_SEVERE_BI_OUT | BO_DEATH | 0.00090000 |
BO_OD_DEATH | BO_OD_DEATH | 1.00000000 |
BO_DEATH | BO_DEATH | 1.00000000 |
BR_ILLICIT | BO_DEATH | 0.01740000 |
BR_ILLICIT | BR_ILLICIT | 999.00000000 |
BR_ILLICIT | BR_OAT_INI | 0.02712532 |
BR_ILLICIT | BR_OD_ILLICIT | 0.03650000 |
BR_OAT_INI | BO_DEATH | 0.01740000 |
BR_OAT_INI | BR_ILLICIT | 999.00000000 |
BR_OAT_INI | BR_OAT_MAINT | 0.84075000 |
BR_OAT_INI | BR_OD_ILLICIT | 0.01285868 |
BR_OAT_MAINT | BO_DEATH | 0.00109940 |
BR_OAT_MAINT | BR_ILLICIT | 999.00000000 |
BR_OAT_MAINT | BR_OAT_MAINT | 0.84075000 |
BR_OAT_MAINT | BR_OD_ILLICIT | 0.01820000 |
BR_OAT_SS | BR_OAT_SS | 999.00000000 |
BR_SS | BR_SS | 999.00000000 |
BR_OD_ILLICIT | BO_MOD_BI | 999.00000000 |
BR_OD_ILLICIT | BO_SEVERE_BI | 0.03000000 |
BR_OD_ILLICIT | BO_OD_DEATH | 0.07700000 |
flextable(costs_tbl) |>
autofit()Variable | Cost |
|---|---|
BN_PN | 0.00 |
BN_ACUTE | 200.34 |
BN_CHRONIC | 70.93 |
BN_CANCER | 580.18 |
BN_OTHER | 70.93 |
BN_PALLIATIVE | 5,405.14 |
BPO_ACUTE | 15,002.06 |
BPO_CHRONIC | 103.05 |
BPO_CANCER | 2,186.73 |
BPO_OTHER | 103.05 |
BPO_PALLIATIVE | 2,904.79 |
BPO_MISUSE | 912.27 |
BI_ILLICIT | 97.10 |
BS_DETOX | 4,550.00 |
BS_OAT_INI | 1,500.88 |
BS_OAT_MAINT | 1,393.88 |
BS_OAT_SS | 1,432.77 |
BS_SS | 1,780.69 |
BO_OD_RX | 974.48 |
BO_OD_ILLICIT | 974.48 |
BO_MOD_BI | 6,825.46 |
BO_SEVERE_BI | 46,749.90 |
BO_SEVERE_BI_OUT | 1,295.96 |
BO_OD_DEATH | 0.00 |
BO_DEATH | 0.00 |
BR_ILLICIT | 97.10 |
BR_OAT_INI | 1,500.88 |
BR_OAT_MAINT | 1,393.88 |
BR_OAT_SS | 1,432.77 |
BR_SS | 1,780.69 |
BR_OD_ILLICIT | 97.10 |
library(here)
library(RColorBrewer)
source(here("02_scripts/01_fun_data.R"))
calib_target_tbl <- read_excel(here("01_data/markov_model_parameters_preprior.xlsx"),
sheet = "calibration_targets")
theme_set(theme_minimal())# prevalence of people on prescribed opioids
v_yrs_prev_rx <- c(2015:2018)
yrs_prev_rx <- length(v_yrs_prev_rx)
prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = yrs_prev_rx,
list(NULL,
paste0("prop_opioids_rx_",
str_sub(v_yrs_prev_rx, 3, 4)))))
tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = yrs_prev_rx,
list(NULL,
paste0("tot_pop_",
str_sub(v_yrs_prev_rx, 3, 4)))))
for (i in 1:yrs_prev_rx) {
yr <- v_yrs_prev_rx[1] + i
prop_opioids_rx_uncalib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_uncalib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(v_yrs_prev_rx,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target, group) %>%
rename(grp = group,
target_val = target) %>%
mutate(year_mon = paste(year, month.abb[7], sep = "_"))
prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Model") %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, target_val) %>%
mutate(grp = "Target"))
prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
levels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
levels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"))
p1 <- ggplot() +
geom_point(data = prop_opioids_rx_uncalib_tbl,
aes(x = year_mon, y = target_val,
color = factor(year), fill = factor(year))) +
geom_point(data = prop_opioids_rx_target_tbl,
aes(x = year_mon, y = target_val),
color = "red") +
geom_point(data = prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon),
aes(x = year_mon, y = target_val),
color = "blue")+
xlab("Year_Month") + ylab("Prevalence of prescription opioid use") +
scale_fill_discrete(name = "year")+
scale_color_discrete(name = "year")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
plotly::ggplotly(p1)The red points represent the observed target proportions of Rx opioids from CIHI report, the blue points represent an annual average (weighted by population during that month), the other coloured points represent monthly prevalence predicted from the model
p2 <- ggplot() +
geom_point(data = prop_opioids_rx_uncalib_wei_mean,
aes(x = year, y = target_val, color = factor(grp),
fill = factor(grp))) +
xlab("Year") + ylab("Prevalence of prescription opioid use")+
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# ylim(0,0.2) +
# ggtitle("Comparison between uncalibrated annual weighted average and
# observed targets of prevalence of Rx opioids use")
plotly::ggplotly(p2)################ Deaths
# number of deaths
v_yr1_deaths <- c(2016:2020)
yrs_deaths <- length(v_yr1_deaths)
num_deaths_uncalib <- rep(NA, yrs_deaths)
for (i in 1:yrs_deaths){
yr <- (v_yr1_deaths[1] - 1) + i
num_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
}
# number of OD-deaths
v_yr1_oddeaths <- c(2016:2021)
yrs_oddeaths <- length(v_yr1_oddeaths)
num_od_deaths_uncalib <- rep(NA, yrs_oddeaths)
for (i in 1:yrs_oddeaths){
yr <- (v_yr1_oddeaths[1] - 1) + i
num_od_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
deaths_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_deaths", "total_od_deaths")) %>%
mutate(target = ifelse(target == "total_deaths", "Total deaths",
ifelse(target == "total_od_deaths",
"Total opioid-related overdose deaths", NA)),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_deaths_uncalib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
mutate(target = "Total opioid-related overdose deaths")) %>%
mutate(year = c(v_yr1_deaths, v_yr1_oddeaths),
group = "Model"))
p3 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl,
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "") +
facet_wrap(~target, scales = "free") +
theme(legend.position="bottom")
plotly::ggplotly(p3)plotly::ggplotly(ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>%
filter(target == "Total deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")) # facet_wrap(~target, scales = "free") + theme(legend.position="bottom")
plotly::ggplotly(ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>%
filter(target == "Total opioid-related overdose deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total opioid-related overdose deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = ""))# number of oat
v_yr1_oat <- c(2018:2021)
yrs_oat <- length(v_yr1_oat)
num_oat_uncalib <- rep(NA, yrs_oat)
for (i in 1:yrs_oat){
yr <- (v_yr1_oat - 1) + i
num_oat_uncalib[i] <- sum(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}
oat_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% "total_oat") %>%
mutate(target = ifelse(target == "total_oat",
"Total OAT", NA),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
mutate(target = "Total OAT",
year = v_yr1_oat,
group = "Model"))
p3 <- ggplot() +
geom_point(data = oat_target_uncalib_tbl,
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "") +
theme(legend.position="bottom")
plotly::ggplotly(p3)# build-in color palette
set.seed(122345)
colors_pal <- Polychrome::createPalette(n_states,
c("#084C61", "#DB504A",
"#E3B505", "#4F6D7A",
"#56A3A6"))
names(colors_pal) <- NULL
trace_tbl_inc <- as_tibble(mod_basecase$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count")
trace_tbl_inc$count[trace_tbl_inc$state == "BO_OD_DEATH" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_OD_DEATH" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_od_deaths
trace_tbl_inc$count[trace_tbl_inc$state == "BO_MOD_BI" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_MOD_BI" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_mod_bi
trace_tbl_inc$count[trace_tbl_inc$state == "BO_SEVERE_BI" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_SEVERE_BI" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_sev_bi
trace_tbl_inc$state <- factor(trace_tbl_inc$state,
levels = v_state_names,
labels = v_state_names)
p4 <- ggplot(trace_tbl_inc, aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl_inc %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal)
plotly::ggplotly(p4)p5 <- ggplot(trace_tbl_inc %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = c(2015:2029)),
aes(x = year, y = count)) +
geom_line() +
ylab("Total opioid-related overdose deaths") + xlab("Year")
# +
# ggtitle("Total opioid-related overdose deaths over time")
plotly::ggplotly(p5)It takes into account fentanyl contamination
library(here)
source(here("02_scripts/03a_fun_calib_td.R"))
theme_set(theme_minimal())
opt_params <- readRDS(here("01_data/map_calib_td.RDS"))Used overall deaths, prevalence of rx opioid use, and total OAT as calibration targets
mod_opt_map <- mod_calib(as.numeric(opt_params$par))# prevalence of people on prescribed opioids
v_yrlast_prev_rx <- c(2015:2018)
yrs_prev_rx <- length(v_yrlast_prev_rx)
prop_opioids_rx_calib <- prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = yrs_prev_rx,
list(NULL,
paste0("prop_opioids_rx_",
str_sub(v_yrlast_prev_rx, 3, 4)))))
tot_pop_calib_tbl <- tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = yrs_prev_rx,
list(NULL,
paste0("tot_pop_",
str_sub(v_yrlast_prev_rx, 3, 4)))))
for (i in 1:yrs_prev_rx) {
yr <- (v_yrlast_prev_rx[1] - 1) + i
prop_opioids_rx_uncalib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_uncalib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
prop_opioids_rx_calib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_calib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(2015:2018,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_calib_tbl <- prop_opioids_rx_calib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(2015:2018,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_calib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target, group) %>%
rename(grp = group,
target_val = target) %>%
mutate(year_mon = paste(year, month.abb[6], sep = "_"))
prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Uncalibrated Model") %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, target_val) %>%
mutate(grp = "target")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Calibrated Model"))
prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_calib_tbl$year_mon <- factor(prop_opioids_rx_calib_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_tbl <- prop_opioids_rx_uncalib_tbl %>%
mutate(group = "Uncalibrated Model") %>%
bind_rows(.,prop_opioids_rx_calib_tbl %>%
mutate(group = "Calibrated Model"))
wprop_opioids_rx_tbl <- prop_opioids_rx_target_tbl %>%
mutate(grp = "Target") %>%
bind_rows(., prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
grp = "Uncalibrated Model")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
grp = "Calibrated Model"))
prop_opioids_rx_tbl$group <- factor(prop_opioids_rx_tbl$group,
labels = c("Calibrated Model", "Uncalibrated Model"),
levels = c("Calibrated Model", "Uncalibrated Model"))
p1 <- ggplot() +
geom_point(data = prop_opioids_rx_tbl, aes(x = year_mon, y = target_val,
color = group, fill = group)) +
geom_point(data = wprop_opioids_rx_tbl %>% filter(grp == "Target"),
aes(x = year_mon, y = target_val, color = grp, fill = grp), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Year_Month") + ylab("Prevalence of Rx opioid use") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
plotly::ggplotly(p1)wprop_opioids_rx_tbl$grp <- factor(wprop_opioids_rx_tbl$grp,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p2 <- ggplot() +
geom_point(data = wprop_opioids_rx_tbl,
aes(x = year_mon, y = target_val,
color = grp, fill = grp), size = 1.5) +
xlab("Year") + ylab("Prevalence of Rx opioid use")+
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# ylim(0,0.2) +
plotly::ggplotly(p2)################ Deaths
# number of deaths
num_deaths_calib <- num_deaths_uncalib <- rep(NA, length(2016:2020))
for (i in 1:length(2016:2020)){
yr <- 2015 + i
num_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
num_deaths_calib[i] <- mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
}
# number of OD-deaths
num_od_deaths_calib <- num_od_deaths_uncalib <- rep(NA, length(2016:2021))
for (i in 1:length(2016:2021)){
yr <- 2015 + i
num_od_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
num_od_deaths_calib[i] <- mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
deaths_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_deaths", "total_od_deaths")) %>%
mutate(target = ifelse(target == "total_deaths", "Total deaths",
ifelse(target == "total_od_deaths",
"Total OD-deaths", NA)),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_deaths_uncalib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
mutate(target = "Total OD-deaths")) %>%
mutate(year = c(2016:2020, 2016:2021),
group = "Uncalibrated Model")) %>%
bind_rows(., data.frame(target_val = num_deaths_calib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_calib) %>%
mutate(target = "Total OD-deaths")) %>%
mutate(year = c(2016:2020, 2016:2021),
group = "Calibrated Model"))
deaths_target_uncalib_tbl$group <- factor(deaths_target_uncalib_tbl$group,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p3 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# +
# facet_wrap(~target, scales = "free") + theme(legend.position="bottom")
plotly::ggplotly(p3)p4 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total OD-deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Opioid-related Overdose Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# +
# facet_wrap(~target, scales = "free") + theme(legend.position="bottom")
plotly::ggplotly(p4)# number of oat
num_oat_calib <- num_oat_uncalib <- rep(NA, length(2018:2021))
for (i in 1:length(2018:2021)){
yr <- 2017 + i
num_oat_uncalib[i] <- sum(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
num_oat_calib[i] <- sum(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}
oat_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_oat")) %>%
mutate(target = ifelse(target == "total_oat",
"Total OAT", NA),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
mutate(target = "Total OAT",
year = c(2018:2021),
group = "Uncalibrated Model")) %>%
bind_rows(., data.frame(target_val = num_oat_calib) %>%
mutate(target = "Total OAT",
year = c(2018:2021),
group = "Calibrated Model"))
oat_target_uncalib_tbl$group <- factor(oat_target_uncalib_tbl$group,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p5 <- ggplot() +
geom_point(data = oat_target_uncalib_tbl,
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total OAT counts") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# +
# facet_wrap(~target, scales = "free") + theme(legend.position="bottom")
plotly::ggplotly(p5)cbind.data.frame(parameter = calib_param_tbl$var_params,
basevalue = calib_param_tbl$basevalue,
opt_map_nm = opt_params$par) %>%
filter(basevalue != opt_map_nm) parameter basevalue opt_map_nm
1 p__BN_PN__BPO_MISUSE 0.00032511 0.00220011
2 p__BPO_CHRONIC__BO_OD_RX 0.00097898 0.00160398
3 p__BPO_PALLIATIVE__BO_DEATH 0.37500000 0.38250000
4 p__BS_DETOX__BI_ILLICIT 0.18181118 0.22181118
5 p__BS_DETOX__BS_OAT_INI 0.80000000 0.82000000
6 p__BS_OAT_MAINT__BS_OAT_MAINT 0.70000000 0.71000000
It doesn’t take into account fentanyl contamination
library(here)
source(here("02_scripts/03b_fun_calib_ntd.R"))
theme_set(theme_minimal())
opt_params_notimedepend <- readRDS(here("01_data/map_calib_ntd.RDS"))Used overall deaths, prevalence of rx opioid use, and total OAT as calibration targets
mod_opt_map_ntd <- mod_calib_ntd(as.numeric(opt_params_notimedepend$par))# prevalence of people on prescribed opioids
prop_opioids_rx_calib <- prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = length(2015:2018),
list(NULL,
paste0("prop_opioids_rx_",
str_sub(2015:2018, 3, 4)))))
tot_pop_calib_tbl <- tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = length(2015:2018),
list(NULL,
paste0("tot_pop_",
str_sub(2015:2018, 3, 4)))))
for (i in 1:length(2015:2018)) {
yr <- 2014 + i
prop_opioids_rx_uncalib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_uncalib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
prop_opioids_rx_calib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_calib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(2015:2018,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_calib_tbl <- prop_opioids_rx_calib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(2015:2018,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_calib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target, group) %>%
rename(grp = group,
target_val = target) %>%
mutate(year_mon = paste(year, month.abb[6], sep = "_"))
prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Uncalibrated Model") %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, target_val) %>%
mutate(grp = "target")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Calibrated Model"))
prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_calib_tbl$year_mon <- factor(prop_opioids_rx_calib_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_tbl <- prop_opioids_rx_uncalib_tbl %>%
mutate(group = "Uncalibrated Model") %>%
bind_rows(.,prop_opioids_rx_calib_tbl %>%
mutate(group = "Calibrated Model"))
wprop_opioids_rx_tbl <- prop_opioids_rx_target_tbl %>%
mutate(grp = "Target") %>%
bind_rows(., prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
grp = "Uncalibrated Model")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
grp = "Calibrated Model"))
prop_opioids_rx_tbl$group <- factor(prop_opioids_rx_tbl$group,
labels = c("Calibrated Model", "Uncalibrated Model"),
levels = c("Calibrated Model", "Uncalibrated Model"))
p1 <- ggplot() +
geom_point(data = prop_opioids_rx_tbl, aes(x = year_mon, y = target_val,
color = group, fill = group)) +
geom_point(data = wprop_opioids_rx_tbl %>% filter(grp == "Target"),
aes(x = year_mon, y = target_val, color = grp, fill = grp), size = 2) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Year_Month") + ylab("Prevalence of Rx opioid use") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
plotly::ggplotly(p1)wprop_opioids_rx_tbl$grp <- factor(wprop_opioids_rx_tbl$grp,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p2 <- ggplot() +
geom_point(data = wprop_opioids_rx_tbl,
aes(x = year_mon, y = target_val,
color = grp, fill = factor(grp)), size = 2) +
xlab("Year") + ylab("Prevalence of Rx opioid use") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# ylim(0,0.2) +
plotly::ggplotly(p2)################ Deaths
# number of deaths
num_deaths_calib <- num_deaths_uncalib <- rep(NA, length(2016:2020))
for (i in 1:length(2016:2020)){
yr <- 2015 + i
num_deaths_uncalib[i] <- mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
num_deaths_calib[i] <- mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
}
# number of OD-deaths
num_od_deaths_calib <- num_od_deaths_uncalib <- rep(NA, length(2016:2021))
for (i in 1:length(2016:2021)){
yr <- 2015 + i
num_od_deaths_uncalib[i] <- mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
num_od_deaths_calib[i] <- mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
deaths_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_deaths", "total_od_deaths")) %>%
mutate(target = ifelse(target == "total_deaths", "Total deaths",
ifelse(target == "total_od_deaths",
"Total OD-deaths", NA)),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_deaths_uncalib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
mutate(target = "Total OD-deaths")) %>%
mutate(year = c(2016:2020, 2016:2021),
group = "Uncalibrated Model")) %>%
bind_rows(., data.frame(target_val = num_deaths_calib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_calib) %>%
mutate(target = "Total OD-deaths")) %>%
mutate(year = c(2016:2020, 2016:2021),
group = "Calibrated Model"))
deaths_target_uncalib_tbl$group <- factor(deaths_target_uncalib_tbl$group,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p3 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# +
# facet_wrap(~target, scales = "free") + theme(legend.position="bottom")
plotly::ggplotly(p3)p4 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total OD-deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Opioid-related Overdose Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
plotly::ggplotly(p4)# number of oat
num_oat_calib <- num_oat_uncalib <- rep(NA, length(2018:2021))
for (i in 1:length(2018:2021)){
yr <- 2017 + i
num_oat_uncalib[i] <- sum(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
num_oat_calib[i] <- sum(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}
oat_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_oat")) %>%
mutate(target = ifelse(target == "total_oat",
"Total OAT", NA),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
mutate(target = "Total OAT",
year = c(2018:2021),
group = "Uncalibrated Model")) %>%
bind_rows(., data.frame(target_val = num_oat_calib) %>%
mutate(target = "Total OAT",
year = c(2018:2021),
group = "Calibrated Model"))
oat_target_uncalib_tbl$group <- factor(oat_target_uncalib_tbl$group,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p5 <- ggplot() +
geom_point(data = oat_target_uncalib_tbl,
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total OAT counts") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# +
# facet_wrap(~target, scales = "free") + theme(legend.position="bottom")
plotly::ggplotly(p5)cbind.data.frame(parameter = calib_param_tbl_ntd$var_params,
basevalue = calib_param_tbl_ntd$basevalue,
opt_map = opt_params_notimedepend$par) %>%
filter(basevalue != opt_map) parameter basevalue opt_map
1 p__BN_PN__BPO_MISUSE 0.00032511 0.00220011
2 p__BPO_CHRONIC__BO_OD_RX 0.00097898 0.00160398
3 p__BPO_PALLIATIVE__BO_DEATH 0.37500000 0.38250000
4 p__BS_OAT_MAINT__BS_OAT_MAINT 0.70000000 0.77000000
library(here)
library(plotly)
source(here("02_scripts/01_fun_data.R"))
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("03_outputs/03_map_calibmod_vs_targets_td.R"))
calib_samples <- readRDS(here("01_data/calib_samples_sir.RDS"))
uncalib_samples <- readRDS(here("01_data/uncalib_sample_sir.RDS"))
opt_params_sir_prev <- readRDS(here("01_data/prev_outcome_data_sir.RDS"))
opt_params_sir_deaths <- readRDS(here("01_data/death_outcome_data_sir.RDS"))
opt_params_sir_oddeaths <- readRDS(here("01_data/oddeath_outcome_data_sir.RDS"))
opt_params_sir_oat <- readRDS(here("01_data/oat_outcome_data_sir.RDS"))
params_sir_prev_uncalib <- readRDS(here("01_data/prev_outcome_data_sir_uncalib.RDS"))
params_sir_deaths_uncalib <- readRDS(here("01_data/death_outcome_data_sir_uncalib.RDS"))
params_sir_oddeaths_uncalib <- readRDS(here("01_data/oddeath_outcome_data_sir_uncalib.RDS"))
params_sir_oat_uncalib <- readRDS(here("01_data/oat_sir_uncalib.RDS"))SIR calibration (calib_mod): first i sampled from prior 10000 sample sets, then i used them to calculated outcomes and likelihoods for each sample set, then i resampled from those 10000 with replacement with likelihood weights 10000 samples (unique sets = 6217), and ran the model on all of them and generated distributions for prevalence of rx opioid use, deaths, od deaths, total oat —- likelihood included deaths (gamma distribution), prevalence of rx opioid use (normal distribution instead of binomial), and total OAT
dt_prev <- as.data.frame(params_sir_prev_uncalib) %>%
# rename(yr15 = "2015",
# yr16 = "2016",
# yr17 = "2017",
# yr18 = "2018") %>%
pivot_longer(1:4, names_to = "prev_year", values_to = "value") %>%
mutate(group = "Uncalibrated Sample") %>%
bind_rows(., as.data.frame(opt_params_sir_prev) %>%
# rename(yr15 = "2015",
# yr16 = "2016",
# yr17 = "2017",
# yr18 = "2018") %>%
pivot_longer(1:4, names_to = "prev_year", values_to = "value") %>%
mutate(group = "SIR Calibrated Sample"))
prev_year = c("yr15", "yr16", "yr17", "yr18")
wprop_opioids_rx_tbl <- prop_opioids_rx_target_tbl %>%
rename(value = target_val) %>%
mutate(grp = "Target",
prev_year = ifelse(grepl("15", year_mon),2015,
ifelse(grepl("16", year_mon), 2016,
ifelse(grepl("17", year_mon), 2017,
ifelse(grepl("18", year_mon), 2018, year_mon)))),
prev_year = factor(prev_year)) %>%
select(-year_mon) %>%
bind_rows(., prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(value = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(prev_year = factor(year),
grp = "Uncalibrated Model - Point estimate")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(value = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(prev_year = factor(year),
grp = "MAP Calibrated Model")) %>%
bind_rows(., dt_prev %>%
group_by(prev_year) %>%
summarize(value = mean(value)) %>%
mutate(grp = "SIR Calibrated Model - Mean",
prev_year = factor(prev_year))) %>%
bind_rows(., dt_prev %>%
group_by(prev_year) %>%
summarize(value = median(value)) %>%
mutate(grp = "SIR Calibrated Model - Median",
prev_year = factor(prev_year)))
ggplotly(ggplot(data = dt_prev,
aes(y = value, x = prev_year)) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = wprop_opioids_rx_tbl,
aes(y = value, x = prev_year, color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group))ggplotly(ggplot(data = dt_prev,
aes(x = value)) +
geom_histogram(aes(fill = group), color="#e9ecef",
alpha=0.3, position = 'identity') +
geom_vline(data = wprop_opioids_rx_tbl, aes(xintercept = value, color = grp))+
theme(axis.text.x = element_text(angle = 45)) +
xlab("Prevalence of rx opioids use") + ylab("") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~prev_year, scales = "free"))p_prev_sir <- ggplot(data = dt_prev,
aes(y = value, x = prev_year)) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = wprop_opioids_rx_tbl,
aes(y = value, x = prev_year, color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group)
ggsave(here("04_report/p_prev_sir.png"))dt_deaths <- as.data.frame(params_sir_deaths_uncalib) %>%
pivot_longer(1:5, names_to = "death_year", values_to = "value") %>%
mutate(group = "Uncalibrated Sample") %>%
bind_rows(., as.data.frame(opt_params_sir_deaths) %>%
pivot_longer(1:5, names_to = "death_year", values_to = "value") %>%
mutate(group = "SIR Calibrated Sample"))
ovrall_deaths_target_uncalib_tbl <- deaths_target_uncalib_tbl %>%
filter(target == "Total deaths") %>%
rename(value = target_val) %>%
mutate(grp = ifelse(group == "Calibrated Model", "MAP Calibrated Model",
ifelse(group == "Uncalibrated Model", "Uncalibrated Model - Point estimate",
ifelse(group == "Target", "Target", group))),
death_year = factor(year)) %>%
select(-c(group, year)) %>%
bind_rows(., dt_deaths %>%
group_by(death_year) %>%
summarize(value = mean(value)) %>%
mutate(grp = "SIR Calibrated Model - Mean",
death_year = factor(death_year))) %>%
bind_rows(., dt_deaths %>%
group_by(death_year) %>%
summarize(value = median(value)) %>%
mutate(grp = "SIR Calibrated Model - Median",
death_year = factor(death_year)))
ggplotly(ggplot(data = dt_deaths,
aes(y = value, x = factor(death_year))) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = ovrall_deaths_target_uncalib_tbl,
aes(y = value, x = factor(death_year), color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group))ggplotly(ggplot(data = dt_deaths,
aes(x = value/1000)) +
geom_histogram(aes(fill = group), color="#e9ecef",
alpha=0.4, position = 'identity') +
# scale_fill_manual(values=c("#69b3a2", "#404080")) +
geom_vline(data = ovrall_deaths_target_uncalib_tbl,
aes(xintercept = value/1000, color = grp))+
theme(axis.text.x = element_text(angle = 45)) +
xlab("Total deaths (in thousands)")+ ylab("") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
xlim(NA, 320)+
facet_wrap(~death_year, scales = "free"))p_deaths_sir <- ggplot(data = dt_deaths,
aes(y = value, x = factor(death_year))) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = ovrall_deaths_target_uncalib_tbl,
aes(y = value, x = factor(death_year), color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group)
ggsave(here("04_report/p_deaths_sir.png"))dt_oat <- as.data.frame(params_sir_oat_uncalib) %>%
pivot_longer(1:length(2018:2021), names_to = "oat_year", values_to = "value") %>%
mutate(group = "Uncalibrated Sample") %>%
bind_rows(., as.data.frame(opt_params_sir_oat) %>%
pivot_longer(1:length(2018:2021), names_to = "oat_year", values_to = "value") %>%
mutate(group = "SIR Calibrated Sample"))
oat_target_uncalib_tbl <- oat_target_uncalib_tbl %>%
rename(value = target_val) %>%
mutate(grp = ifelse(group == "Calibrated Model", "MAP Calibrated Model",
ifelse(group == "Uncalibrated Model", "Uncalibrated Model - Point estimate",
ifelse(group == "Target", "Target", group))),
oat_year = factor(year)) %>%
select(-c(group, year)) %>%
bind_rows(., dt_oat %>%
group_by(oat_year) %>%
summarize(value = mean(value)) %>%
mutate(grp = "SIR Calibrated Model - Mean",
oat_year = factor(oat_year))) %>%
bind_rows(., dt_oat %>%
group_by(oat_year) %>%
summarize(value = median(value)) %>%
mutate(grp = "SIR Calibrated Model - Median",
oat_year = factor(oat_year)))
ggplotly(ggplot(data = dt_oat,
aes(y = value, x = factor(oat_year))) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = oat_target_uncalib_tbl,
aes(y = value, x = factor(oat_year), color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group))ggplotly(ggplot(data = dt_oat,
aes(x = value)) +
geom_histogram(aes(fill = group), color="#e9ecef",
alpha=0.4, position = 'identity') +
geom_vline(data = oat_target_uncalib_tbl,
aes(xintercept = value, color = grp))+
theme(axis.text.x = element_text(angle = 45)) +
xlab("Total OAT")+ ylab("") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~oat_year, scales = "free"))p_oat_sir <- ggplot(data = dt_oat,
aes(y = value, x = factor(oat_year))) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = oat_target_uncalib_tbl,
aes(y = value, x = factor(oat_year), color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group)
ggsave(here("04_report/p_oat_sir.png"))library(here)
library(flextable)
library(plotly)
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("02_scripts/06_interventions_models.R"))
theme_set(theme_minimal())Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 Scenario 2: Scenario 2: Doesn’t account for Fentanyl contamination and has the same probability from Illicit use to OD and from OD to OD death
trace_tbl <- as_tibble(mod_basecase$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "No Interventions") %>%
bind_rows(., as_tibble(mod_nalx_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Naloxone")) %>%
bind_rows(., as_tibble(mod_ss_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Safer Supply")) %>%
bind_rows(., as_tibble(mod_pres_guid_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Prescription Guidelines")) %>%
bind_rows(., as_tibble(mod_all_int_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "All Interventions"))
trace_tbl$mod <- factor(trace_tbl$mod,
levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))
v_extra_sevbi <- c(mod_basecase$extra_sev_bi, mod_nalx_1$extra_sev_bi,
mod_ss_1$extra_sev_bi, mod_pres_guid_1$extra_sev_bi,
mod_all_int_1$extra_sev_bi)
v_extra_modbi <- c(mod_basecase$extra_mod_bi, mod_nalx_1$extra_mod_bi,
mod_ss_1$extra_mod_bi, mod_pres_guid_1$extra_mod_bi,
mod_all_int_1$extra_mod_bi)
v_extra_od <- c(mod_basecase$extra_od_deaths, mod_nalx_1$extra_od_deaths,
mod_ss_1$extra_od_deaths, mod_pres_guid_1$extra_od_deaths,
mod_all_int_1$extra_od_deaths)
mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)
for (j in 1:length(mod_names)){
trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_od[j]
trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_modbi[j]
trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_sevbi[j]
}
trace_tbl$state <- factor(trace_tbl$state,
levels = v_state_names,
labels = v_state_names)trace_tbl_2 <- as_tibble(mod_basecase_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "No Interventions") %>%
bind_rows(., as_tibble(mod_nalx_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Naloxone")) %>%
bind_rows(., as_tibble(mod_ss_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Safer Supply")) %>%
bind_rows(., as_tibble(mod_pres_guid_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Prescription Guidelines")) %>%
bind_rows(., as_tibble(mod_all_int_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "All Interventions"))
trace_tbl_2$mod <- factor(trace_tbl_2$mod,
levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))
v_extra_sevbi_2 <- c(mod_basecase_2$extra_sev_bi, mod_nalx_2$extra_sev_bi,
mod_ss_2$extra_sev_bi, mod_pres_guid_2$extra_sev_bi,
mod_all_int_2$extra_sev_bi)
v_extra_modbi_2 <- c(mod_basecase_2$extra_mod_bi, mod_nalx_2$extra_mod_bi,
mod_ss_2$extra_mod_bi, mod_pres_guid_2$extra_mod_bi,
mod_all_int_2$extra_mod_bi)
v_extra_od_2 <- c(mod_basecase_2$extra_od_deaths, mod_nalx_2$extra_od_deaths,
mod_ss_2$extra_od_deaths, mod_pres_guid_2$extra_od_deaths,
mod_all_int_2$extra_od_deaths)
for (j in 1:length(mod_names)){
trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_od_2[j]
trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_modbi_2[j]
trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_sevbi_2[j]
}
trace_tbl_2$state <- factor(trace_tbl_2$state,
levels = v_state_names,
labels = v_state_names)
trace_tbl <- trace_tbl %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., trace_tbl_2 %>%
mutate(scenario = "Scenario 2"))p1 <- ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 1") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p1)p2 <- ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 2") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p2)plotly::ggplotly(ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = rep(c(2015:2029), length(mod_names)*2)),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention") +
facet_wrap(~scenario) +
labs(caption = "Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 \n
Scenario 2: Doesn't account for Fentanyl contamination and has the same probability from \n Illicit use to OD and from OD to OD death"))inci_od_deaths_basecase <- inci_od_deaths_nalox <- inci_od_deaths_ss <- inci_od_deaths_pg <- inci_od_deaths_all <- inci_od_deaths_basecase_2 <- inci_od_deaths_nalox_2 <- inci_od_deaths_ss_2 <- inci_od_deaths_pg_2 <- inci_od_deaths_all_2 <- rep(NA, 14)
for (i in 1:14){
yr <- 2015 + i
inci_od_deaths_basecase[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_basecase_2[i] <- mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_nalox[i] <- mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_nalox_2[i] <- mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_ss[i] <- mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_ss_2[i] <- mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_pg[i] <- mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_pg_2[i] <- mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_all[i] <- mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_all_2[i] <- mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
inc_od_death_tbl <- data.frame(intervention = c(rep(mod_names, each = 14)),
year = c(rep(2016:2029, 5)),
inci_od_deaths = c(inci_od_deaths_basecase,
inci_od_deaths_nalox,
inci_od_deaths_ss,
inci_od_deaths_pg,
inci_od_deaths_all),
scenario = rep("Scenario 1", 14*5)) %>%
bind_rows(., data.frame(intervention = c(rep(mod_names, each = 14)),
year = c(rep(2016:2029, 5)),
inci_od_deaths = c(inci_od_deaths_basecase_2,
inci_od_deaths_nalox_2,
inci_od_deaths_ss_2,
inci_od_deaths_pg_2,
inci_od_deaths_all_2),
scenario = rep("Scenario 2", 14*5)))
saveRDS(inc_od_death_tbl, file = here("01_data/inc_od_death_tbl_mod.RDS"))
target_od_deaths <- calib_target_tbl %>%
filter(group == "total_od_deaths") %>%
select(year, target) %>%
rename(inci_od_deaths = target) %>%
mutate(intervention = "Target")
saveRDS(target_od_deaths, file = here("01_data/inc_od_death_tbl_target.RDS"))
inc_od_death_tbl$intervention <- factor(inc_od_death_tbl$intervention,
levels = mod_names, labels = mod_names)p3 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 1"),
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2")
ggplotly(p3)p4 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 2"),
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2")
ggplotly(p4)plotly::ggplotly(ggplot(inc_od_death_tbl,
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2") +
facet_wrap(~scenario))p5 <- ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 1") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p5)p6 <- ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 2") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p6)plotly::ggplotly(ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = rep(c(2015:2029), length(mod_names)*2)),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention") +
facet_wrap(~scenario))set.seed(122345)
colors_pal <- Polychrome::createPalette(n_states,
c("#084C61", "#DB504A",
"#E3B505", "#4F6D7A",
"#56A3A6"))
names(colors_pal) <- NULL
ggplotly(ggplot(trace_tbl %>%
filter(grepl("BPO", state)), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>% filter(cycle_num == 180) %>%
filter(grepl("BPO", state)),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BPO_MISUSE", "BI_ILLICIT",
"BS_OAT_INI", "BS_OAT_MAINT")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BPO_MISUSE", "BI_ILLICIT",
"BS_OAT_INI", "BS_OAT_MAINT")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BS_OAT_SS", "BS_SS")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BS_OAT_SS", "BS_SS")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
"BO_MOD_BI", "BO_SEVERE_BI")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
"BO_MOD_BI", "BO_SEVERE_BI")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
"BR_OAT_MAINT", "BR_OAT_SS",
"BR_SS", "BR_OD_ILLICIT")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
"BR_OAT_MAINT", "BR_OAT_SS",
"BR_SS", "BR_OD_ILLICIT")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))cost_diff <- c(mod_nalx_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_ss_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_pres_guid_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_all_int_1$total_net_present_cost - mod_basecase$total_net_present_cost)
cost_diff_per <- c(round((cost_diff/c(mod_basecase$total_net_present_cost))*100,2))
deaths_eff_tbl <- trace_tbl %>%
filter(scenario == "Scenario 1") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>%
mutate(nalx_effect = Naloxone - `No Interventions`,
ss_effect = `Safer Supply` - `No Interventions`,
pg_effect = `Prescription Guidelines` - `No Interventions`,
all_effect = `All Interventions` - `No Interventions`,
nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>%
select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))
total_death_oddeaths <- trace_tbl %>%
filter(scenario == "Scenario 1") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH"))
od_deaths_diff <- c(deaths_eff_tbl$nalx_effect[1], deaths_eff_tbl$ss_effect[1],
deaths_eff_tbl$pg_effect[1], deaths_eff_tbl$all_effect[1])
od_deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[1], deaths_eff_tbl$ss_effect_per[1],
deaths_eff_tbl$pg_effect_per[1], deaths_eff_tbl$all_effect_per[1])
deaths_diff <- c(deaths_eff_tbl$nalx_effect[2], deaths_eff_tbl$ss_effect[2],
deaths_eff_tbl$pg_effect[2], deaths_eff_tbl$all_effect[2])
deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[2], deaths_eff_tbl$ss_effect_per[2],
deaths_eff_tbl$pg_effect_per[2], deaths_eff_tbl$all_effect_per[2])
saveRDS(cbind.data.frame(cost_diff, cost_diff_per,
od_deaths_diff, od_deaths_diff_per,
deaths_diff, deaths_diff_per), here("01_data/point_estiamte.RDS"))
saveRDS(total_death_oddeaths, here("01_data/total_deaths_oddeaths.RDS"))
costs <- paste0(round(cost_diff/1000000, 0), " (", cost_diff_per, "%)")
deaths <- paste0(round(deaths_diff, 0), " (", deaths_diff_per, "%)")
oddeaths <- paste0(round(od_deaths_diff, 0), " (", od_deaths_diff_per, "%)")
effects_tbl <- cbind.data.frame(mod_names[-1], costs, deaths, oddeaths)
effects_tbl |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Mean Change Compared to No Intervention")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
costs = "Discounted Net Present Costs \n in Millions (%)",
deaths = "Deaths (%)",
oddeaths = "Overdose-related Deaths (%)")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Mean Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths (%) | Overdose-related Deaths (%) |
Naloxone | 102 (0.01%) | 396 (0.01%) | -6420 (-4.05%) |
Safer Supply | 4233 (0.53%) | -13956 (-0.31%) | -3138 (-1.98%) |
Prescription Guidelines | -26103 (-3.28%) | 15252 (0.34%) | -1764 (-1.11%) |
All Interventions | -24876 (-3.13%) | 14115 (0.31%) | -8543 (-5.39%) |
cost_diff_2 <- c(mod_nalx_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_ss_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_pres_guid_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_all_int_2$total_net_present_cost - mod_basecase_2$total_net_present_cost)
cost_diff_per_2 <- c(round((cost_diff_2/c(mod_basecase_2$total_net_present_cost))*100,2))
deaths_eff_tbl_2 <- trace_tbl %>%
filter(scenario == "Scenario 2") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>%
mutate(nalx_effect = Naloxone - `No Interventions`,
ss_effect = `Safer Supply` - `No Interventions`,
pg_effect = `Prescription Guidelines` - `No Interventions`,
all_effect = `All Interventions` - `No Interventions`,
nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>%
select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))
od_deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[1], deaths_eff_tbl_2$ss_effect[1],
deaths_eff_tbl_2$pg_effect[1], deaths_eff_tbl_2$all_effect[1])
od_deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[1], deaths_eff_tbl_2$ss_effect_per[1],
deaths_eff_tbl_2$pg_effect_per[1], deaths_eff_tbl_2$all_effect_per[1])
deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[2], deaths_eff_tbl_2$ss_effect[2],
deaths_eff_tbl_2$pg_effect[2], deaths_eff_tbl_2$all_effect[2])
deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[2], deaths_eff_tbl_2$ss_effect_per[2],
deaths_eff_tbl_2$pg_effect_per[2], deaths_eff_tbl_2$all_effect_per[2])
costs_2 <- paste0(round(cost_diff_2/1000000, 0), " (", cost_diff_per_2, "%)")
deaths_2 <- paste0(round(deaths_diff_2, 0), " (", deaths_diff_per_2, "%)")
oddeaths_2 <- paste0(round(od_deaths_diff_2, 0), " (", od_deaths_diff_per_2, "%)")
effects_tbl_2 <- cbind.data.frame(mod_names[-1], costs_2, deaths_2, oddeaths_2)
effects_tbl_2 |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Mean Change Compared to No Intervention")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
costs_2 = "Discounted Net Present Costs \n in Millions (%)",
deaths_2 = "Deaths (%)",
oddeaths_2 = "Overdose-related Deaths (%)")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Mean Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths (%) | Overdose-related Deaths (%) |
Naloxone | 85 (0.01%) | 303 (0.01%) | -4861 (-3.9%) |
Safer Supply | 4047 (0.51%) | -12743 (-0.28%) | -2435 (-1.95%) |
Prescription Guidelines | -26098 (-3.29%) | 15308 (0.34%) | -1485 (-1.19%) |
All Interventions | -24882 (-3.13%) | 14105 (0.31%) | -6607 (-5.3%) |
tbl_effect_plot1 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
deaths_diff_per, od_deaths_diff_per) %>%
pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
deaths_diff_per_2, od_deaths_diff_per_2) %>%
pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>%
mutate(scenario = "Scenario 2")) %>%
mutate(group = ifelse(group == "cost_diff_per_2", "cost_diff_per",
ifelse(group == "deaths_diff_per_2", "deaths_diff_per",
ifelse(group == "od_deaths_diff_per_2", "od_deaths_diff_per", group))))p7 <- ggplot(data = tbl_effect_plot1 %>%
filter(scenario == "Scenario 1"), aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_minimal()
ggplotly(p7)p8 <- ggplot(data = tbl_effect_plot1 %>%
filter(scenario == "Scenario 2"), aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_minimal()
ggplotly(p8)ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
# scale_shape_manual(values = c(15, 16, 17, 18)) +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_bw() +
facet_wrap(~scenario))tbl_effect_plot2 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
od_deaths_diff_per) %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
od_deaths_diff_per_2) %>%
rename(cost_diff_per = "cost_diff_per_2",
od_deaths_diff_per = "od_deaths_diff_per_2") %>%
mutate(scenario = "Scenario 2"))p9 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 1"),
aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_minimal()
ggplotly(p9)p10 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 2"),
aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_minimal()
ggplotly(p10)ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_bw() +
facet_wrap(~scenario))v_yrs_prev <- c(2015:2029)
yrs_prev <- length(v_yrs_prev)
prev_fun <- function(mod, i, v_states){
prop_uncalib <- data.frame(matrix(NA,byrow = T, nrow = 12,
ncol = yrs_prev,
list(NULL,
paste0("prop_",
str_sub(v_yrs_prev, 3, 4)))))
tot_pop_uncalib_tbl <- data.frame(matrix(NA,byrow = T, nrow = 12,
ncol = yrs_prev,
list(NULL,
paste0("tot_pop_",
str_sub(v_yrs_prev, 3, 4)))))
for (j in 1:yrs_prev) {
yr <- (v_yrs_prev[1] - 1) + j
if (length(v_states) > 1)
{
prop_uncalib[, j] <- assign(
paste0("prop_", str_sub(yr, 3, 4)),
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
v_states]) /
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
} else {
prop_uncalib[, j] <- assign(
paste0("prop_", str_sub(yr, 3, 4)),
mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
v_states]) /
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
}
tot_pop_uncalib_tbl[, j] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
# final table for monthly prevalence over 15 years
prop_uncalib_tbl <- prop_uncalib %>%
pivot_longer(1:15, names_to = "grp", values_to = "value") %>%
arrange(grp) %>%
mutate(year = rep(v_yrs_prev,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:15, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop)) %>%
mutate(intervention = mod_names[i])
# table with weighted yearly prevalence
prop_uncalib_wei_mean <- prop_uncalib_tbl %>%
group_by(year) %>%
summarise(value = weighted.mean(value, tot_pop_val)) %>%
ungroup() %>%
mutate(intervention = mod_names[i])
return (list(prop_uncalib_tbl = prop_uncalib_tbl,
prop_uncalib_wei_mean = prop_uncalib_wei_mean))
}noint_prev_sevbi_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
noint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
nalx_prev_sevbi_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
nalx_prev_sevbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
ss_prev_sevbi_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
ss_prev_sevbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
pg_prev_sevbi_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
pg_prev_sevbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
allint_prev_sevbi_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
allint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
prev_sevbi_mon_tbl <- noint_prev_sevbi_mon_tbl %>%
bind_rows(., nalx_prev_sevbi_mon_tbl, ss_prev_sevbi_mon_tbl,
pg_prev_sevbi_mon_tbl, allint_prev_sevbi_mon_tbl)
prev_sevbi_yr_tbl <- noint_prev_sevbi_yr_tbl %>%
bind_rows(., nalx_prev_sevbi_yr_tbl, ss_prev_sevbi_yr_tbl,
pg_prev_sevbi_yr_tbl, allint_prev_sevbi_yr_tbl)
prev_sevbi_mon_tbl$year_mon <- factor(prev_sevbi_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BO_SEVERE_BI") %>%
filter(scenario == "Scenario 1") %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))ggplotly(ggplot() +
geom_point(data = prev_sevbi_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of severe brain injury") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))noint_prev_modbi_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
noint_prev_modbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
nalx_prev_modbi_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
nalx_prev_modbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
ss_prev_modbi_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
ss_prev_modbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
pg_prev_modbi_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
pg_prev_modbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
allint_prev_modbi_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
allint_prev_modbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
prev_modbi_mon_tbl <- noint_prev_modbi_mon_tbl %>%
bind_rows(., nalx_prev_modbi_mon_tbl, ss_prev_modbi_mon_tbl,
pg_prev_modbi_mon_tbl, allint_prev_modbi_mon_tbl)
prev_modbi_yr_tbl <- noint_prev_modbi_yr_tbl %>%
bind_rows(., nalx_prev_modbi_yr_tbl, ss_prev_modbi_yr_tbl,
pg_prev_modbi_yr_tbl, allint_prev_modbi_yr_tbl)
prev_modbi_mon_tbl$year_mon <- factor(prev_modbi_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BO_MOD_BI") %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))ggplotly(ggplot() +
geom_point(data = prev_modbi_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of moderate brain injury") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))v_states_prev_tx <- c("BS_DETOX", "BS_OAT_INI", "BS_OAT_MAINT",
"BS_OAT_SS", "BS_SS", "BR_OAT_INI", "BR_OAT_MAINT",
"BR_OAT_SS", "BR_SS")
noint_prev_tx_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_tx)$prop_uncalib_tbl
noint_prev_tx_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
nalx_prev_tx_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_prev_tx)$prop_uncalib_tbl
nalx_prev_tx_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
ss_prev_tx_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_prev_tx)$prop_uncalib_tbl
ss_prev_tx_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
pg_prev_tx_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_prev_tx)$prop_uncalib_tbl
pg_prev_tx_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
allint_prev_tx_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_prev_tx)$prop_uncalib_tbl
allint_prev_tx_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
prev_tx_mon_tbl <- noint_prev_tx_mon_tbl %>%
bind_rows(., nalx_prev_tx_mon_tbl, ss_prev_tx_mon_tbl,
pg_prev_tx_mon_tbl, allint_prev_tx_mon_tbl)
prev_tx_yr_tbl <- noint_prev_tx_yr_tbl %>%
bind_rows(., nalx_prev_tx_yr_tbl, ss_prev_tx_yr_tbl,
pg_prev_tx_yr_tbl, allint_prev_tx_yr_tbl)
prev_tx_mon_tbl$year_mon <- factor(prev_tx_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_prev_tx) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))ggplotly(ggplot() +
geom_point(data = prev_tx_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of substance use treatment") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target) %>%
rename(value = target) %>%
mutate(year_mon = paste(year, month.abb[7], sep = "_"))
v_states_prev_rx_opd <- c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")
noint_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
noint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
nalx_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
nalx_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
ss_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
ss_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
pg_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
pg_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
allint_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
allint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
prev_rx_opd_mon_tbl <- noint_prev_rx_opd_mon_tbl %>%
bind_rows(., nalx_prev_rx_opd_mon_tbl, ss_prev_rx_opd_mon_tbl,
pg_prev_rx_opd_mon_tbl, allint_prev_rx_opd_mon_tbl) %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, value, year_mon) %>%
mutate(intervention = "Target"))
prev_rx_opd_yr_tbl <- noint_prev_rx_opd_yr_tbl %>%
bind_rows(., nalx_prev_rx_opd_yr_tbl, ss_prev_rx_opd_yr_tbl,
pg_prev_rx_opd_yr_tbl, allint_prev_rx_opd_yr_tbl) %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, value) %>%
mutate(intervention = "Target"))
prev_rx_opd_mon_tbl$year_mon <- factor(prev_rx_opd_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
levels = paste(rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))ggplotly(ggplot() +
geom_point(data = prev_rx_opd_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of prescription opioid use") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))v_states_illicit <- c("BI_ILLICIT", "BR_ILLICIT")
noint_prev_illicituse_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicit)$prop_uncalib_tbl
noint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicit)$prop_uncalib_wei_mean
nalx_prev_illicituse_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_illicit)$prop_uncalib_tbl
nalx_prev_illicituse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_illicit)$prop_uncalib_wei_mean
ss_prev_illicituse_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_illicit)$prop_uncalib_tbl
ss_prev_illicituse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_illicit)$prop_uncalib_wei_mean
pg_prev_illicituse_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_illicit)$prop_uncalib_tbl
pg_prev_illicituse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_illicit)$prop_uncalib_wei_mean
allint_prev_illicituse_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_illicit)$prop_uncalib_tbl
allint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_illicit)$prop_uncalib_wei_mean
prev_illicituse_mon_tbl <- noint_prev_illicituse_mon_tbl %>%
bind_rows(., nalx_prev_illicituse_mon_tbl, ss_prev_illicituse_mon_tbl,
pg_prev_illicituse_mon_tbl, allint_prev_illicituse_mon_tbl)
prev_illicituse_yr_tbl <- noint_prev_illicituse_yr_tbl %>%
bind_rows(., nalx_prev_illicituse_yr_tbl, ss_prev_illicituse_yr_tbl,
pg_prev_illicituse_yr_tbl, allint_prev_illicituse_yr_tbl)
prev_illicituse_mon_tbl$year_mon <- factor(prev_illicituse_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_illicit) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))ggplotly(ggplot() +
geom_point(data = prev_illicituse_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of illicit opioid use") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))noint_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
noint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
nalx_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
nalx_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
ss_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
ss_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
pg_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
pg_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
allint_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
allint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
prev_rxmisuse_mon_tbl <- noint_prev_rxmisuse_mon_tbl %>%
bind_rows(., nalx_prev_rxmisuse_mon_tbl, ss_prev_rxmisuse_mon_tbl,
pg_prev_rxmisuse_mon_tbl, allint_prev_rxmisuse_mon_tbl)
prev_rxmisuse_yr_tbl <- noint_prev_rxmisuse_yr_tbl %>%
bind_rows(., nalx_prev_rxmisuse_yr_tbl, ss_prev_rxmisuse_yr_tbl,
pg_prev_rxmisuse_yr_tbl, allint_prev_rxmisuse_yr_tbl)
prev_rxmisuse_mon_tbl$year_mon <- factor(prev_rxmisuse_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BPO_MISUSE") %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))ggplotly(ggplot() +
geom_point(data = prev_rxmisuse_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of rx opioid misuse") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))v_states_od <- c("BO_OD_RX", "BO_OD_ILLICIT", "BR_OD_ILLICIT")
noint_prev_od_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_od)$prop_uncalib_tbl
noint_prev_od_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_od)$prop_uncalib_wei_mean
nalx_prev_od_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_od)$prop_uncalib_tbl
nalx_prev_od_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_od)$prop_uncalib_wei_mean
ss_prev_od_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_od)$prop_uncalib_tbl
ss_prev_od_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_od)$prop_uncalib_wei_mean
pg_prev_od_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_od)$prop_uncalib_tbl
pg_prev_od_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_od)$prop_uncalib_wei_mean
allint_prev_od_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_od)$prop_uncalib_tbl
allint_prev_od_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_od)$prop_uncalib_wei_mean
prev_od_mon_tbl <- noint_prev_od_mon_tbl %>%
bind_rows(., nalx_prev_od_mon_tbl, ss_prev_od_mon_tbl,
pg_prev_od_mon_tbl, allint_prev_od_mon_tbl)
prev_od_yr_tbl <- noint_prev_od_yr_tbl %>%
bind_rows(., nalx_prev_od_yr_tbl, ss_prev_od_yr_tbl,
pg_prev_od_yr_tbl, allint_prev_od_yr_tbl)
prev_od_mon_tbl$year_mon <- factor(prev_od_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_od) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))ggplotly(ggplot() +
geom_point(data = prev_od_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of opioid-related overdose") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))v_states_illicitod <- c("BO_OD_ILLICIT", "BR_OD_ILLICIT")
noint_prev_illicitod_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicitod)$prop_uncalib_tbl
noint_prev_illicitod_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
nalx_prev_illicitod_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_illicitod)$prop_uncalib_tbl
nalx_prev_illicitod_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
ss_prev_illicitod_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_illicitod)$prop_uncalib_tbl
ss_prev_illicitod_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
pg_prev_illicitod_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_illicitod)$prop_uncalib_tbl
pg_prev_illicitod_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
allint_prev_illicitod_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_illicitod)$prop_uncalib_tbl
allint_prev_illicitod_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
prev_illicitod_mon_tbl <- noint_prev_illicitod_mon_tbl %>%
bind_rows(., nalx_prev_illicitod_mon_tbl, ss_prev_illicitod_mon_tbl,
pg_prev_illicitod_mon_tbl, allint_prev_illicitod_mon_tbl)
prev_illicitod_yr_tbl <- noint_prev_illicitod_yr_tbl %>%
bind_rows(., nalx_prev_illicitod_yr_tbl, ss_prev_illicitod_yr_tbl,
pg_prev_illicitod_yr_tbl, allint_prev_illicitod_yr_tbl)
prev_illicitod_mon_tbl$year_mon <- factor(prev_illicitod_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_illicitod) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))ggplotly(ggplot() +
geom_point(data = prev_illicitod_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of illicit opioid-related overdose") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))noint_prev_rxod_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_OD_RX")$prop_uncalib_tbl
noint_prev_rxod_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
nalx_prev_rxod_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BO_OD_RX")$prop_uncalib_tbl
nalx_prev_rxod_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
ss_prev_rxod_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BO_OD_RX")$prop_uncalib_tbl
ss_prev_rxod_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
pg_prev_rxod_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BO_OD_RX")$prop_uncalib_tbl
pg_prev_rxod_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
allint_prev_rxod_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BO_OD_RX")$prop_uncalib_tbl
allint_prev_rxod_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
prev_rxod_mon_tbl <- noint_prev_rxod_mon_tbl %>%
bind_rows(., nalx_prev_rxod_mon_tbl, ss_prev_rxod_mon_tbl,
pg_prev_rxod_mon_tbl, allint_prev_rxod_mon_tbl)
prev_rxod_yr_tbl <- noint_prev_rxod_yr_tbl %>%
bind_rows(., nalx_prev_rxod_yr_tbl, ss_prev_rxod_yr_tbl,
pg_prev_rxod_yr_tbl, allint_prev_rxod_yr_tbl)
prev_rxod_mon_tbl$year_mon <- factor(prev_rxod_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BO_OD_RX") %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))ggplotly(ggplot() +
geom_point(data = prev_rxod_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of rx opioid-related overdose") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))library(tidyverse)
library(plotly)
library(flextable)
theme_set(theme_minimal())
library(here)
m_outcomes_psa <- readRDS(here("01_data/m_outcomes_psa.RDS"))
point_est <- readRDS(here("01_data/point_estiamte.RDS"))
total_death_oddeaths <- readRDS(here("01_data/total_deaths_oddeaths.RDS"))
mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)ggplotly(ggplot(data = diff_psa_dt %>%
filter(grp == "death"), aes(x = diff, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Change in Deaths Compared to No Intervention") + ylab("") +
geom_vline(data = point_est_diff %>%
filter(grp == "death"), aes(xintercept = diff, linetype = type))+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal())ggplotly(ggplot(data = diff_psa_dt %>%
filter(grp == "cost"), aes(x = diff, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Change in Costs Compared to No Intervention (in millions)") + ylab("") +
geom_vline(data = point_est_diff %>%
filter(grp == "cost"), aes(xintercept = diff/1000000, linetype = type), linewidth = 0.5)+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal())diff_per_psa_dt <- outcomes_psa_dt %>%
mutate(cost_diff_per_nalox = cost_diff_per_nalox,
cost_diff_per_ss = cost_diff_per_ss,
cost_diff_per_pg = cost_diff_per_pg,
cost_diff_per_all = cost_diff_per_all) %>%
select(cost_diff_per_nalox, cost_diff_per_ss, cost_diff_per_pg, cost_diff_per_all,
death_diff_per_nalox, death_diff_per_ss, death_diff_per_pg, death_diff_per_all,
oddeath_diff_per_nalox, oddeath_diff_per_ss, oddeath_diff_per_pg, oddeath_diff_per_all) %>%
pivot_longer(1:12, names_to = "group", values_to = "diff_per") %>%
separate(col = group, into = c("grp", "type1", "type2", "Intervention"), sep = "_") %>%
select(-type1, -type2) %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions", Intervention)))))
point_est_diff_per <- point_est %>%
select(cost_diff_per, deaths_diff_per, od_deaths_diff_per) %>%
cbind.data.frame(Intervention = mod_names[-1]) %>%
pivot_longer(1:3,names_to = "grp", values_to = "diff") %>%
separate(grp, into = c("grp", "type"), sep = "_") %>% select(-type) %>%
mutate(grp = ifelse(grp == "od", "oddeath",
ifelse(grp == "deaths", "death", grp)),
type = "Point Estimate")ggplotly(ggplot(data = diff_per_psa_dt %>%
filter(grp == "death"), aes(x = diff_per, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("% Change in Deaths Compared to No Intervention") + ylab("") +
geom_vline(data = point_est_diff_per %>%
filter(grp == "death"), aes(xintercept = diff, linetype = type))+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal())ggplotly(ggplot(data = diff_per_psa_dt %>%
filter(grp == "cost"), aes(x = diff_per, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("% Change in Costs Compared to No Intervention") + ylab("") +
geom_vline(data = point_est_diff_per %>%
filter(grp == "cost"), aes(xintercept = diff, linetype = type))+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal())mean_cost_diff_nalox <- mean(m_outcomes_psa[, "cost_diff_nalox"])
ci_cost_diff_nalox <- quantile(m_outcomes_psa[, "cost_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_ss <- mean(m_outcomes_psa[, "cost_diff_ss"])
ci_cost_diff_ss <- quantile(m_outcomes_psa[, "cost_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_pg <- mean(m_outcomes_psa[, "cost_diff_pg"])
ci_cost_diff_pg <- quantile(m_outcomes_psa[, "cost_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_all <- mean(m_outcomes_psa[, "cost_diff_all"])
ci_cost_diff_all <- quantile(m_outcomes_psa[, "cost_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff <- c(paste0(round(ci_cost_diff_nalox[3]/1000000, 0), " [",
round(ci_cost_diff_nalox[1]/1000000, 0),", ",
round(ci_cost_diff_nalox[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_ss[3]/1000000, 0), " [",
round(ci_cost_diff_ss[1]/1000000, 0),", ",
round(ci_cost_diff_ss[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_pg[3]/1000000, 0), " [",
round(ci_cost_diff_pg[1]/1000000, 0),", ",
round(ci_cost_diff_pg[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_all[3]/1000000, 0), " [",
round(ci_cost_diff_all[1]/1000000, 0),", ",
round(ci_cost_diff_all[2]/1000000, 0), "]"))
mean_death_diff_nalox <- mean(m_outcomes_psa[, "death_diff_nalox"])
ci_death_diff_nalox <- quantile(m_outcomes_psa[, "death_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_ss <- mean(m_outcomes_psa[, "death_diff_ss"])
ci_death_diff_ss <- quantile(m_outcomes_psa[, "death_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_pg <- mean(m_outcomes_psa[, "death_diff_pg"])
ci_death_diff_pg <- quantile(m_outcomes_psa[, "death_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_all <- mean(m_outcomes_psa[, "death_diff_all"])
ci_death_diff_all <- quantile(m_outcomes_psa[, "death_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff <- c(paste0(round(ci_death_diff_nalox[3], 0), " [",
round(ci_death_diff_nalox[1], 0),", ",
round(ci_death_diff_nalox[2], 0), "]"),
paste0(round(ci_death_diff_ss[3], 0), " [",
round(ci_death_diff_ss[1], 0),", ",
round(ci_death_diff_ss[2], 0), "]"),
paste0(round(ci_death_diff_pg[3], 0), " [",
round(ci_death_diff_pg[1], 0),", ",
round(ci_death_diff_pg[2], 0), "]"),
paste0(round(ci_death_diff_all[3], 0), " [",
round(ci_death_diff_all[1], 0),", ",
round(ci_death_diff_all[2], 0), "]"))
mean_oddeath_diff_nalox <- mean(m_outcomes_psa[, "oddeath_diff_nalox"])
ci_oddeath_diff_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_ss <- mean(m_outcomes_psa[, "oddeath_diff_ss"])
ci_oddeath_diff_ss <- quantile(m_outcomes_psa[, "oddeath_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_pg <- mean(m_outcomes_psa[, "oddeath_diff_pg"])
ci_oddeath_diff_pg <- quantile(m_outcomes_psa[, "oddeath_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_all <- mean(m_outcomes_psa[, "oddeath_diff_all"])
ci_oddeath_diff_all <- quantile(m_outcomes_psa[, "oddeath_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff <- c(paste0(round(ci_oddeath_diff_nalox[3], 0), " [",
round(ci_oddeath_diff_nalox[1], 0),", ",
round(ci_oddeath_diff_nalox[2], 0), "]"),
paste0(round(ci_oddeath_diff_ss[3], 0), " [",
round(ci_oddeath_diff_ss[1], 0),", ",
round(ci_oddeath_diff_ss[2], 0), "]"),
paste0(round(ci_oddeath_diff_pg[3], 0), " [",
round(ci_oddeath_diff_pg[1], 0),", ",
round(ci_oddeath_diff_pg[2], 0), "]"),
paste0(round(ci_oddeath_diff_all[3], 0), " [",
round(ci_oddeath_diff_all[1], 0),", ",
round(ci_oddeath_diff_all[2], 0), "]"))
effects_tbl <- cbind.data.frame(mod_names[-1], mean_cost_diff,
mean_death_diff, mean_oddeath_diff)
effects_tbl |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
mean_cost_diff = "Discounted Net Present \n Costs in Millions ",
mean_death_diff = "Deaths ",
mean_oddeath_diff = "Opioid-related Overdose Deaths")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present | Deaths | Opioid-related Overdose Deaths |
Naloxone | 209 [107, 362] | 625 [414, 926] | -10089 [-15009, -6697] |
Safer Supply | 4216 [3651, 4821] | -14180 [-17532, -11074] | -3221 [-4105, -2529] |
Prescription Guidelines | -25455 [-31480, -19866] | 14417 [8610, 20565] | -5458 [-11045, -1624] |
All Interventions | -21058 [-27059, -15403] | 960 [-5687, 7753] | -18544 [-28496, -11309] |
mean_cost_diff_per_nalox <- mean(m_outcomes_psa[, "cost_diff_per_nalox"])
ci_cost_diff_per_nalox <- quantile(m_outcomes_psa[, "cost_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_ss <- mean(m_outcomes_psa[, "cost_diff_per_ss"])
ci_cost_diff_per_ss <- quantile(m_outcomes_psa[, "cost_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_pg <- mean(m_outcomes_psa[, "cost_diff_per_pg"])
ci_cost_diff_per_pg <- quantile(m_outcomes_psa[, "cost_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_all <- mean(m_outcomes_psa[, "cost_diff_per_all"])
ci_cost_diff_per_all <- quantile(m_outcomes_psa[, "cost_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per <- c(paste0(round(ci_cost_diff_per_nalox[3], 2), " [",
round(ci_cost_diff_per_nalox[1], 2),", ",
round(ci_cost_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_cost_diff_per_ss[3], 2), " [",
round(ci_cost_diff_per_ss[1], 2),", ",
round(ci_cost_diff_per_ss[2], 2), "]%"),
paste0(round(ci_cost_diff_per_pg[3], 2), " [",
round(ci_cost_diff_per_pg[1], 2),", ",
round(ci_cost_diff_per_pg[2], 2), "]%"),
paste0(round(ci_cost_diff_per_all[3], 2), " [",
round(ci_cost_diff_per_all[1], 2),", ",
round(ci_cost_diff_per_all[2], 2), "]%"))
mean_death_diff_per_nalox <- mean(m_outcomes_psa[, "death_diff_per_nalox"])
ci_death_diff_per_nalox <- quantile(m_outcomes_psa[, "death_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_ss <- mean(m_outcomes_psa[, "death_diff_per_ss"])
ci_death_diff_per_ss <- quantile(m_outcomes_psa[, "death_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_pg <- mean(m_outcomes_psa[, "death_diff_per_pg"])
ci_death_diff_per_pg <- quantile(m_outcomes_psa[, "death_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_all <- mean(m_outcomes_psa[, "death_diff_per_all"])
ci_death_diff_per_all <- quantile(m_outcomes_psa[, "death_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per <- c(paste0(round(ci_death_diff_per_nalox[3], 2), " [",
round(ci_death_diff_per_nalox[1], 2),", ",
round(ci_death_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_death_diff_per_ss[3], 2), " [",
round(ci_death_diff_per_ss[1], 2),", ",
round(ci_death_diff_per_ss[2], 2), "]%"),
paste0(round(ci_death_diff_per_pg[3], 2), " [",
round(ci_death_diff_per_pg[1], 2),", ",
round(ci_death_diff_per_pg[2], 2), "]%"),
paste0(round(ci_death_diff_per_all[3], 2), " [",
round(ci_death_diff_per_all[1], 2),", ",
round(ci_death_diff_per_all[2], 2), "]%"))
mean_oddeath_diff_per_nalox <- mean(m_outcomes_psa[, "oddeath_diff_per_nalox"])
ci_oddeath_diff_per_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_ss <- mean(m_outcomes_psa[, "oddeath_diff_per_ss"])
ci_oddeath_diff_per_ss <- quantile(m_outcomes_psa[, "oddeath_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_pg <- mean(m_outcomes_psa[, "oddeath_diff_per_pg"])
ci_oddeath_diff_per_pg <- quantile(m_outcomes_psa[, "oddeath_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_all <- mean(m_outcomes_psa[, "oddeath_diff_per_all"])
ci_oddeath_diff_per_all <- quantile(m_outcomes_psa[, "oddeath_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per <- c(paste0(round(ci_oddeath_diff_per_nalox[3], 2), " [",
round(ci_oddeath_diff_per_nalox[1], 2),", ",
round(ci_oddeath_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_ss[3], 2), " [",
round(ci_oddeath_diff_per_ss[1], 2),", ",
round(ci_oddeath_diff_per_ss[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_pg[3], 2), " [",
round(ci_oddeath_diff_per_pg[1], 2),", ",
round(ci_oddeath_diff_per_pg[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_all[3], 2), " [",
round(ci_oddeath_diff_per_all[1], 2),", ",
round(ci_oddeath_diff_per_all[2], 2), "]%"))
effects_tbl_per <- cbind.data.frame(mod_names[-1], mean_cost_diff_per,
mean_death_diff_per, mean_oddeath_diff_per)
effects_tbl_per |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Percent Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
mean_cost_diff_per = "Discounted Net Present Costs",
mean_death_diff_per = "Deaths ",
mean_oddeath_diff_per = "Opioid-related Overdose Deaths")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Percent Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths | Opioid-related Overdose Deaths |
Naloxone | 0.03 [0.01, 0.05]% | 0.01 [0.01, 0.02]% | -3.46 [-4.13, -3.13]% |
Safer Supply | 0.54 [0.46, 0.61]% | -0.31 [-0.38, -0.24]% | -1.11 [-1.97, -0.66]% |
Prescription Guidelines | -3.23 [-3.94, -2.55]% | 0.32 [0.19, 0.45]% | -1.81 [-2.76, -0.93]% |
All Interventions | -2.67 [-3.39, -1.98]% | 0.02 [-0.13, 0.17]% | -6.43 [-7.16, -5.58]% |
tbl_effect_plot1 <- cbind.data.frame(value = c("lb", "ub", "median"),
ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
ci_cost_diff_per_pg, ci_cost_diff_per_all,
ci_death_diff_per_nalox, ci_death_diff_per_ss,
ci_death_diff_per_pg, ci_death_diff_per_all,
ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>%
pivot_longer(2:13, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3)) %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions", Intervention)))))
tbl_effect_plot1$Intervention <- factor(tbl_effect_plot1$Intervention, levels = mod_names, labels = mod_names)
ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = median))+
geom_jitter(aes(color = Intervention, shape = Intervention), size = 2, position = position_dodge(0.25)) +
# geom_errorbar(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), width = 0.3, position = position_dodge(0.1)) +
geom_pointrange(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), position = position_dodge(0.25)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16, 17, 18)) +
# scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","Opioid-related Overdose Deaths"))+
theme_minimal())tbl_effect_plot2 <- cbind.data.frame(value = c("cost_lb", "cost_ub", "cost_median"),
ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
ci_cost_diff_per_pg, ci_cost_diff_per_all) %>%
pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3, group)) %>%
full_join(., cbind.data.frame(value = c("oddeath_lb", "oddeath_ub", "oddeath_median"),
ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>%
pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3, group)), by = "Intervention") %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions", Intervention)))))
tbl_effect_plot2$Intervention <- factor(tbl_effect_plot2$Intervention, levels = mod_names, labels = mod_names)
ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_median, y = oddeath_median)) +
geom_point(aes(color = Intervention), size = 1) +
geom_pointrange(aes(ymin = oddeath_lb, ymax = oddeath_ub, color = Intervention)) +
geom_errorbarh(aes(xmax = cost_lb, xmin = cost_ub,color = Intervention), height = 0) +
scale_color_brewer(palette = "Dark2") +
# scale_shape_manual(values = c(15, 16, 17, 18)) +
geom_hline(yintercept = 0, linewidth = 0.25) +
geom_vline(xintercept = 0, linewidth = 0.25) +
xlab("Change in Costs Compared to No Intervention (%)") +
ylab("Change in Opioid-related Overdose Deaths \n Compared to No Intervention (%)") +
theme_minimal())death_psa_dt <- outcomes_psa_dt %>%
select(tot_death_no, tot_death_nalox, tot_death_ss, tot_death_pg, tot_death_all,
tot_oddeath_no, tot_oddeath_nalox, tot_oddeath_ss, tot_oddeath_pg, tot_oddeath_all) %>%
pivot_longer(1:10, names_to = "group", values_to = "total") %>%
separate(col = group, into = c("type", "grp", "Intervention"), sep = "_") %>%
select(-type) %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions",
ifelse(Intervention == "no", "No Interventions", Intervention))))))
death_psa_dt$Intervention <- factor(death_psa_dt$Intervention,
levels = mod_names, labels = mod_names)
total_death_oddeaths <- total_death_oddeaths %>%
pivot_longer(3:7, names_to = "Intervention", values_to = "total") %>%
mutate(grp = ifelse(state == "BO_OD_DEATH", "oddeath", "death"),
type = "Point Estimate") %>%
select(-state)
ggplot(data = death_psa_dt %>%
filter(grp == "oddeath"), aes(x = total, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Total Opioid-related Overdose Deaths after 15 years") + ylab("") +
geom_vline(data = total_death_oddeaths %>%
filter(grp == "oddeath"), aes(xintercept = total, linetype = type), linewidth = 0.5)+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal()ggplot(data = death_psa_dt %>%
filter(grp == "death"), aes(x = total/1000, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Total Deaths after 15 years (in thousands)") + ylab("") +
geom_vline(data = total_death_oddeaths %>%
filter(grp == "death"), aes(xintercept = total/1000, linetype = type), linewidth = 0.5)+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal()outcomes_psa_dt_inci <- outcomes_psa_dt %>%
select(35:104)
inci_od_dt <- data.frame(nm = rep(NA, length(outcomes_psa_dt_inci)),
med = rep(NA, length(outcomes_psa_dt_inci)),
lb = rep(NA, length(outcomes_psa_dt_inci)),
ub = rep(NA, length(outcomes_psa_dt_inci)))
for (i in 1:length(outcomes_psa_dt_inci)){
inci_od_dt$nm[i] <- names(outcomes_psa_dt_inci)[i]
inci_od_dt[i, 2:4] <- as.numeric(quantile(outcomes_psa_dt_inci[i, ], c(0.5, 0.025, 0.975)))
}
inc_od_death_tbl_mod <- readRDS(here("01_data/inc_od_death_tbl_mod.RDS")) %>%
filter(scenario == "Scenario 1") %>%
select(-scenario) %>%
rename(Intervention = intervention,
Year = year,
value = inci_od_deaths) %>%
mutate(grp = "Point Estimate")
inci_od_dt_plot <- inci_od_dt %>%
separate(nm, into = c("type", "type1", "Intervention", "Year"), sep = "_") %>%
select(-c(type, type1)) %>%
rename(value = med) %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions",
ifelse(Intervention == "no", "No Interventions", Intervention))))),
Year = c(rep(2016:2029, 5)),
grp = "PSA results") %>%
bind_rows(., inc_od_death_tbl_mod) %>%
mutate(newgrp = paste0(Intervention, "_", grp))
inci_od_dt_plot$Intervention <- factor(inci_od_dt_plot$Intervention, levels = mod_names, labels = mod_names)
inc_od_death_tbl_target <- readRDS(here("01_data/inc_od_death_tbl_target.RDS")) %>%
mutate(grp = "Target")
ggplotly(ggplot(data = inci_od_dt_plot, aes(x = Year, y = value, color = grp, fill = grp)) +
geom_line() +
geom_point(data = inc_od_death_tbl_target, aes(x = year, y = inci_od_deaths, color = grp, fill = grp))+
scale_color_manual(values = RColorBrewer::brewer.pal(3, "Dark2")[c(1,3,2)]) +
scale_fill_manual(values = RColorBrewer::brewer.pal(3, "Dark2")[c(1,3, 2)]) +
# scale_color_brewer(palette = "Dark2") +
# scale_fill_brewer(palette = "Dark2") +
geom_ribbon(aes(ymin = lb, ymax = ub), alpha = 0.5) +
facet_wrap(~Intervention, scales = "free") +
theme(legend.title = element_blank()))